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Introduction 


Broadly  speaking  our  work  aimed  to  develop  and  assess  the  utility  of  various  theoretical 
and  experimental  diffuse  light  imaging  schemes  for  tumor  detection  and  characterization 
within  the  human  breast.  During  the  last  year  we  have  theoretically  and  numerically 
investigated  the  angular  spectrum  approach  in  the  parallel  plate  (compressed  breast) 
geometry.  On  the  experimental  side  we  have  built  a  soft  compression  plate  apparatus  for 
diffuse  optical  tomography  of  the  breast.  The  device  was  used  to  experimentally  to  test 
the  angular  spectrum  and  other  theoretical  approaches.  We  have  pilot  tested  the 
instrument  with  tissue  phantoms  and  in  a  few  clinical  scenarios.  Through  this  work  we 
uncovered  some  technical  limitations  of  diffuse  optical  tomography  based  on  the  angular 
spectrum  approach,  and,  importantly,  we  developed  other  schemes  for  improved  breast 
imaging. 

Body 

Below  we  have  summarized  our  progress  related  to  each  of  three  specific  aims  outlined 
for  this  grant. 

Specific  aim  1 

To  theoretically  explore  the  utility  of  the  FFT  approach  for  breast  tumor  imaging  with 
diffusing  light  in  the  parallel  plate  compressed  breast  configuration. 

During  the  last  year  we  contributed  2  papers  that  developed  diffraction  tomographic 
techniques  for  diffuse  optical  imaging.  The  first  paper  investigated  the  use  of  depth 
dependent  regularization  and  mathematical  filtering  methods  to  improve  the  fidelity  of 
these  projection  images  in  the  presence  of  noise  [1].  In  This  paper  we  developed  an 
algorithm  that  stabilized  the  diffraction  tomographic  deconvolution  for  optimal  resolution 
at  all  depths.  In  addition  we  found  an  empirical  method  to  estimate  the  depth  of  a 
heterogeneity  of  interest  using  only  2D  projection  images.  In  a  related  publication  we 
documented  a  group  of  refinements  to  the  diffraction  tomographic  techniques  [2],  These 
improvements  included  incorporating  boundary  conditions  in  order  to  extend  the 
technique  to  slab  media  of  finite  size  and  to  semi-infinite  media.  The  modified  theory 
with  better  boundary  conditions  was  demonstrated  to  be  superior  to  the  initial  formulation 
using  numerical  simulations  plus  noise.  We  also  studied  the  cross-coupling  of  absorption 
and  scattering;  an  important  finding  of  these  numerical  experiments  was  that  cross¬ 
coupling  of  absorption  and  scattering  would  arise  in  samples  for  which  the  thin  slice 
approximation  breaks  down.  This  finding  suggested  that  2D  projection  images  alone  will 
not  be  sufficient  to  reconstruct  true  optical  properties  in  many  practical  situations,  and 
that  some  fonn  of  3D  reconstruction  is  critical. 

Specific  aim  2 

To  construct  an  apparatus  based  on  photon  migration  technology  that  will  enable  us  to 
test  image  reconstruction  approaches  with  realistic  conditions  and  signal  to  noise  ratios. 
We  have  made  substantial  progress  in  instrumentation  development  during  the  last  year. 
There  were  several  important  issues  to  be  faced  as  we  began  this  project.  Since  the  focus 
of  the  grant  was  on  the  FFT  approach,  we  designed  a  compressed  breast  instrument  for 
which  it  was  possible  to  rigorously  test  these  theoretical  schemes.  At  the  same  time 
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however,  we  wanted  to  design  an  apparatus  that  would  be  adaptable  to  other  theoretical 
and  experimental  approaches  that  we  hoped  to  investigate. 

So  far  we  have  primarily  been  using  the  frequency  domain  version  of  the  apparatus  to 
measure  the  optical  properties  of  tissue  phantoms  and  of  normal  breast  [3],  At  the  heart  of 
the  device  are  two  parallel  compression  plates  (e.g.  a  “source-plate”  and  an  “output- 
plate”)  with  adjustable  separation.  The  original  version  of  the  source  plate  had  one 
source  fiber  coupled  in  series  to  3  amplitude  modulated  (140MHz)  diode  lasers  operating 
at  3  wavelengths  in  the  NIR  (750nm,  786nm,  and  830nm).  The  original  version  of  the 
detection  plate  had  a  scannable  detection  fiber.  Typically  we  would  scan  the  fiber  over 
an  -lOxlOcm  area  in  steps  of  0.5  cm,  which  required  about  20  minutes 

We  used  the  original  device  to  demonstrate  the  basic  ideas  of  the  angular  spectrum 
approach,  and  also  to  show  the  limitations  of  the  approach  experimentally  in  tissue 
phantoms.  We  also  used  the  device  in  pilot  measurements  of  the  average  optical 
properties  of  normal  human  breast.  In  carrying  out  these  latter  studies  we  obtained  a 
better  understanding  of  some  of  the  clinical  problems  that  will  arise  in  practice.  For 
example,  proper  treatment  of  the  chest  boundary  was  critical  for  obtaining  correction 
optical  properties.  It  also  became  apparent  that  one  needed  to  obtain  a  larger  data  sets, 
more  quickly  (i.e.  than  20  minutes)  in  order  to  push  the  methodology  effectively  into  the 
clinic,  and  to  carry  out  realistic  3D  image  reconstructions.  In  this  vein  we  have  designed 
and  started  building  a  second-generation  optical  mammography  device  that  employs  a 
CCD  camera  that  rapidly  collects  large  data  sets  [4], 

Specific  aim  3 

To  theoretically  formulate  and  demonstrate  non-perturbative  approaches  to  homogenous 
and  heterogenous  media  that  explicitly  incorporates  the  boundary  conditions  and  that  is 
accurate  when  variations  in  scattering  and  absorption  are  large  (i.e.  not  in  the 
perturbative  limit). 

In  a  this  vein  we  derived  a  new  integro-differential  equation  for  bounded  media  within 
the  diffusion  approximation  (also  following  the  principles  of  diffraction  tomography) 
which  readily  makes  possible  full  3D  reconstructions  in  the  compressed  breast  geometry 
[5].  The  new  technique  employs  a  series  of  plane  diffuse  photon  density  waves  with 
different  modulation  frequencies.  It  requires  both  two-dimensional  FFT  operations  (as  in 
the  original  approach)  and  a  one-dimensional  matrix  inversion. 

In  addition  we  developed  finite  difference  numerical  codes  for  data  inversion. 
This  approach  was  utilized  to  obtain  measurements  of  the  bulk  optical  properties  of  tissue 
phantoms  and  normal  breast  using  the  apparatus  developed  in  specific  aim  2.  During  the 
course  of  analyzing  the  results  of  clinical  measurements  it  became  apparent  that  our 
analytical  schemes  were  suffering  in  practical  scenarios.  As  a  result  we  developed  finite 
difference  numerical  codes  for  data  inversion.  These  numerical  schemes  can  flexibly 
segment  and  extend  the  regions  of  interest.  In  particular  we  have  been  using  the  explicit 
adjoint  formulation  for  the  inverse  problem.  Most  DOT  approaches  have  been 
implemented  in  2D.  In  the  best  2D  cases,  researchers  have  employed  a  cylindrical 
geometry  to  reduce  the  dimensionality  of  the  problem.  In  fact,  only  recently  have  a  few 
direct  comparisons  between  full  3D  and  2D  reconstructions  been  made  using  simulated 
data;  not  surprisingly  the  3D  reconstructions  were  superior.  Furthermore  in  any  clinical 
setting  the  forward  problem  is  necessarily  3D.  This  is  especially  true  for  our  compressed 
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plate  geometry  where  there  is  no  intrinsic  symmetry  about  which  to  reduce  the 

dimensionality  of  the  problem.  A  description  of  these  numerical  schemes  is  provided  [6]. 

Key  Research  Accomplishments 

•  Developed  an  algorithm  that  stabilized  the  diffraction  tomographic  technique  for 
optimal  resolution  at  all  depths. 

•  Found  an  empirical  method  to  estimate  the  depth  of  a  heterogeneity  of  interest  using 
only  2D  projection  images. 

•  Incorporated  boundary  conditions  into  the  diffraction  tomography  approach  in  order 
to  extend  the  technique  to  slab  media  of  finite  size  and  to  semi-infinite  media. 

•  Derived  a  new  integro-differential  equation  for  bounded  media  within  the  diffusion 
approximation  (also  following  the  principles  of  diffraction  tomography),  which 
readily  makes  possible  full  3D  reconstructions  in  the  compressed  breast  geometry. 

•  Demonstrate  the  basic  ideas  and  tested  the  limitations  of  the  angular  spectrum 
approach  using  clinical  prototype  imager. 

•  Designed  and  started  building  a  second-generation  optical  mammography  device 
employing  a  CCD  camera  to  rapidly  collects  large  data  sets. 

•  Developed  finite  difference  numerical  codes  for  data  inversion. 

Reportable  outcomes 

•  Grant:  Parallel,  Rapid  Diffuse  Optical  Tomography  of  Breast 

•  Grant:  Parallel  Detection  and  Computation  for  Diffuse  Tomography  of  Breast 

•  Durduran,  T.,  Culver,  J.P.,  Hoboke,  M.J.,  Li,  X.D.,  Zubkov,  L.,  Chance,  B., 
Pattanayak,  D.N.,  and  Yodh,  A.G.,  Algorithms  for  3D  localization  and  imaging  using 
near-field  diffraction  tomography  with  diffuse  light,  Optics  Express  4,  247-262 
(1999). 

•  Li,  X.D.,  Pattanayak,  D.N.,  Durduran,  T.,  Culver,  J.P.,  Chance,  B.,  and  Yodh,  A.G., 
Near-field  diffraction  tomography  with  diffuse  photon  density  waves,  (Accepted  by 
Physical  Review  B)  (2000). 

•  Durduran,  T.,  Holboke,  M.,  Culver,  J.P.,  Zubkov,  L.,  Choe,  R.,  Pattanayak,  D.N., 
Chance,  B.,  and  Yodh,  A.G.,  Tissue  bulk  optical  properties  of  breasts  and  phantoms 
obtained  with  clinical  optical  imager  (To  be  published.  Technical  Digest  Biomedical 
Topical  Meeting,  OSA)  (2000). 

•  Culver,  J.P.,  Ntziachristos,  V.,  Zubkov,  L.,  Durduran,  T.,  Pattanayak,  D.N.,  and 
Yodh,  A.G.,  Data  set  size  and  image  quality  in  diffuse  optical  mammography: 
evaluation  of  a  clinical  prototype  (To  be  published,  Technical  Digest  Biomedical 
Topical  Meeting,  OSA)  (2000). 

•  Pattanayak,  D.N.,  and  Yodh,  A.G.,  Diffuse  optical  3D-slice  imaging  of  bounded 
turbid  media  using  a  new  integro-differential  equation.  Optics  Express  4,  231-240 
(1999). 

•  Holboke,  M.J.  and  Yodh,  A.G.,  Parallel  three-dimensional  diffuse  optical  tomography 
(To  be  published,  Technical  Digest  Biomedical  Topical  Meeting,  OSA)  (2000). 

Conclusions 

In  our  work  over  the  last  year  we  developed  and  assessed  the  utility  of  various  theoretical 

and  experimental  diffuse  light  imaging  schemes  for  tumor  detection  and  characterization 
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within  the  human  breast.  Focusing  on  the  parallel  plate  (compressed  breast)  geometry,  we 
have  theoretically  and  numerically  investigated  the  angular  spectrum  approach. 

Significant  progress  was  made  in  stablizing  the  techinque  with  respect  to  measurement 
noise  and  extending  the  technique  to  clinically  relevant  geometries.  Through  this  work 
we  uncovered  some  technical  limitations  of  diffuse  optical  tomography  based  on  the 
angular  spectrum  approach,  and,  importantly,  we  developed  other  schemes  for  improved 
breast  imaging.  In  particular  we  have  derived  a  new  integro-differential  equation  for 
bounded  media  within  the  diffusion  approximation  and  developed  finite  difference 
numerical  codes  for  arbitrary  geometries. 

On  the  experimental  side  we  have  built  a  soft  compression  plate  apparatus  for  diffuse 
optical  tomography  of  the  breast  and  conducted  pilot  tests  of  the  instrument  with  tissue 
phantoms  and  in  a  few  clinical  scenarios.  While  using  the  current  system  for  bulk  tissue 
properties  it  became  obvious  that  the  data  acquisition  speed  needed  to  be  increase 
significantly.  As  a  result  we  designed  and  began  building  a  second-generation 
mammography  machine  using  a  CCD  camera.  The  larger  data  sets  will  be  complemented 
by  future  efforts  to  develop  parallel  processor  image  reconstruction  code. 

In  conclusion  we  have  gained  significant  knowledge  and  experience  with  the  optical 
mammography  problem  during  this  period.  As  a  result  of  this  understanding,  we  have 
designed  and  started  the  development  of  novel  computational  and  experimental  tools  that 
permit  rapid  acquisition  and  analysis  of  informationally  dense  diffuse  optical  data  sets  in 
the  compressed  breast  geometry. 
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Abstract:  We  introduce  two  filtering  methods  for  near-field  diffuse 
light  diffraction  tomography  based  on  the  angular  spectrum  represen¬ 
tation.  We  then  combine  these  filtering  techniques  with  a  new  method 
to  find  the  approximate  depth  of  the  image  heterogeneities.  Taken  to¬ 
gether  these  ideas  improve  the  fidelity  of  our  projection  image  recon¬ 
structions,  provide  an  interesting  three  dimensional  rendering  of  the 
reconstructed  volume,  and  enable  us  to  identify  and  classify  image  ar¬ 
tifacts  that  need  to  be  controlled  better  for  tissue  applications.  The 
analysis  is  accomplished  using  data  derived  from  numerical  finite  dif¬ 
ference  simulations  with  added  noise. 
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1.  Introduction 


Diffuse  photon  density  waves  (DPDWs)  [1-7]  and  their  time-domain  analogs  [8-10] 
can  provide  quantitative  spectroscopic  information  about  tissue  structure,  tissue  chro- 
mophores,  and  tissue  metabolism.  Tomographic  imaging  of  deep  tissues  with  diffusive 
waves  offers  the  possibility  for  functional  imaging  based  on  these  spectroscopic  con¬ 
trasts.  Thus  image  reconstruction  has  been  explored  intensely  with  varying  degrees  of 
success  in  tissue  phantoms  [11,12],  the  brain  [13-18],  and  the  breast  [19-24].  One  theo¬ 
retical  approach  that  has  received  significant  recent  interest  is  the  method  of  diffraction 
tomography  [25-28]  extended  to  the  near- field  [29-43].  This  scheme  is  particularly  well 
suited  for  planar  source-detection  geometries  that  arise  for  example,  in  breast  imaging. 
The  methodology  has  attracted  attention  [29-43]  because  the  diffusive  wave  modes  de¬ 
rived  using  the  angular-spectrum  representation  provide  a  very  convenient  framework 
for  imaging  in  parallel-plate  geometries,  and  for  analytic  or  semi-analytic  investigation 
of  resolution  and  signal-to-noise  issues.  The  method  also  offers  a  very  fast  and  direct 
way  to  obtain  two-dimensional  projection  images  of  the  turbid  media;  these  images  are 
essentially  optical  mammograms  [21]. 

In  this  contribution  we  study  a  set  of  theoretical  issues  associated  with  image 
reconstruction  using  near-field  diffraction  tomography.  In  order  to  obtain  quality  im¬ 
ages,  for  example,  one  apply  spectral  filters  [41-46]  to  the  data  at  several  levels  of  the 
image  processing  (i.e.  filters  with  respect  to  spatial  frequencies  in  the  reconstructions). 
Well  defined  rules  do  not  exist  for  choosing  these  filter  functions;  in  fact  different  filter 
functions  introduce  uncontrolled  systematic  errors  into  the  images  [41]  .  Additionally, 
the  images  are  dramatically  improved  when  the  experimenter  has  knowledge  of  the  ap¬ 
proximate  depth  of  the  optical  heterogeneity.  Object  depth  determination  is  generally 
difficult  unless  one  has  other  means  to  obtain  this  information  such  as  multiple  optical 
views  of  the  same  sample  or  image  segmentation  from  other  types  of  diagnostics. 

We  investigate  both  critical  issues  in  this  paper,  and  demonstrate  algorithms 
for  optimizing  projection  image  fidelity.  Data  for  these  studies  is  derived  from  numerical 
simulations  with  added  noise.  Two  kinds  of  mathematical  filters  are  introduced;  (1)  a 
phase-only  filter  which  does  not  have  any  free  parameters  and  allows  accurate  localiza¬ 
tion,  and  (2)  a  depth  and  noise  dependent  low-pass  filter  with  the  cut-off  frequency  as 
a  free  parameter.  Then  the  two  filtering  methods  are  combined  with  a  technique  to  find 
the  approximate  depth  of  the  heterogeneities.  We  obtain  two  dimensional  projection 
images  and  we  demonstrate  a  three  dimensional  rendering  of  these  image  projections 
which  appears  promising  for  clinical  applications.  Our  work  clarifies  the  mentioned  is¬ 
sues  in  a  quantitative  way;  it  makes  apparent  the  limitations  of  the  technique,  identifies 
artifacts  and  directions  for  improvement. 

2.  Theory 

2. 1  Angular  spectrum  representation  of  diffuse  photon  density  waves 

The  angular-spectrum  representation  provides  a  set  of  modes  well  suited  for  the  descrip¬ 
tion  of  the  propagation  and  scattering  of  diffuse  photon  density  waves  in  parallel-plate 
geometries  (fig.(l)).  The  representation  and  its  applications  have  been  described  by 
numerous  investigators  in  optics  [25-28,47-48]  and  in  the  photon  migration  community 
[29-43].  For  clarity  we  review  the  salient  features  of  this  approach  below.  We  assume 
infinite  space  boundary  conditions  in  this  discussion;  Green’s  functions  for  semi-infinite 
and  slab  geometries  have  been  reported  in  the  literature  [33-35]  and  do  not  qualitatively 
change  the  conclusions  presented  herein.  The  starting  point  for  our  treatment  is  the  for¬ 
mal  expression  relating  the  scattered  diffuse  light  field,  $sc  ,  to  a  volume  integral  over 
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heterogeneities  within  the  entire  turbid  medium  [33]: 

$sc(r)=  [  T(r')G0(r,r')d3r' . 
Jv 


z„ 


(Xci,yd,Zd) 


Detector 
Scans 
in  a  plane 


Figure  1.  The  generic  near-field  diffraction  tomography  experiment.  The  detector 
is  scanned  in  a  2D  grid  on  the  surface  of  the  plane  parallel  to  and  displaced  from 
the  plane  containing  the  source.  The  breast  is  embedded  in  the  box  along  with 
Intralipid  in  order  to  match  its  average  optical  properties. 


Here  T(r')  is  “tumor”  function  at  position  r'  which,  in  the  Born  Approximation, 
depends  on  the  background  diffuse  photon  density  wave,  4>o(r),  as  well  as  the  absorption 
(5^a(r'))  and  scattering  (<5/4( r')  )  deviations  from  the  medium’s  background  absorption 
(pa0  )  and  scattering  coefficients  (/i'0  )  respectively.  T(r')  also  depends  on  the  speed 
of  light  in  the  medium,  u,  the  background  light  diffusion  constant  Do,  the  background 
diffuse  photon  density  wavevector  ko  =  y/(—vfiao  +  ito)/Do)  ,  and  the  modulation  fre¬ 
quency,  c j.  It  can  be  conveniently  separated  into  absorptive  and  scattering  parts  [49], 
i.e., 

Tabs{r )  =  -^-$o(r )<5Ma(r)  ,  (2) 

Do 

Tsc( r)  =  ^Uo(r)  •  (3) 

Here  we  drop  the  gradient  terms  due  to  the  scattering  variations  for  simplicity.  This 
approximation  is  described  in  detail  in  [49].  Go(r,r')  is  the  usual  Green’s  function  for 
the  medium  of  the  form: 

<4) 

Typical  experiments  that  employ  the  angular-spectrum  representation  for  anal¬ 
ysis  use  a  “parallel-plate”  geometry  (e.g.  fig.(l))  where  the  detected  DPDW,  4>(r)  — 
<hsc(r)  +  $0(r)  ,  is  measured  at  discrete  points  within  a  square  grid  on  the  planar  output 
(detector  plane,  (x,y,  za))  surface.  We  take  our  source  (at  the  origin,  (0,0,0))  to  be  a 
point  emitter  on  the  planar  input  surface  (source  plane,  (xs,ys,  0)  ).  In  this  geometry 
the  Green’s  function  is  conveniently  written  in  the  following  form: 

Go(r,  r')  =  J  y’^G0(p,5!2d,^)e“i(p(x""x')+<7(!/,J“s/'))  ,  (5) 

where  the  z-axis  is  defined  in  the  direction  normal  to  the  plane  surfaces,  and  (p,  q )  are 
the  2D  spatial  (k-space)  frequencies  conjugate  to  the  x-  and  ^-coordinates.  The  angular 
spectrum  (Weyl  expansion)  representation  of  the  Green’s  function  is  [31,50] 

Go(p,q,z<i,z')  —  t  (6) 
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where  m  =  y/k^  -  {p2  +  g^)  is  the  complex  propagator  wavenumber  and  the  imaginary 
part  of  m  is  positive,  i.e.  S(m)  >  0.  With  these  definitions,  the  2D  Fourier  Transform 
of  the  scattered  field  in  the  detection  plane  is  simply  [33] 


$sc(p,  q,  Zd )  =  J  dz'Go(p,  q,  zd,  z')T{p,  q,  z)  . 


(7) 


where  hat  (A)  indicates  quantities  in  k-space.  This  equation  is  the  fundamental  equation 
for  diffraction  tomography. 

2.2  Image  reconstruction 

Equation  (7)  relates  the  tumor  function  in  each  slice  to  the  2D  Fourier  Transform  of  the 
scattered  DPDW  in  the  detection  plane.  Rather  than  inverting  the  Laplace  Transform 
directly,  we  discretize  the  integral  into  a  sum  over  slices  of  thickness  Az,  i.e 


N  N  -A 

*,cfaq,zd)  =  Y,Go(p^zd,z')t{p^Zj)  =  £ \JLf(p,q,Zjym^-^  .  (8) 


3  =  1 


J=1 


A z  is  the  step  size  in  the  z-direction  and  N  =  zd/Az  is  the  total  number  of  slices.  This 
is  the  key  equation  for  our  image  reconstruction  using  the  angular-spectrum  approach. 
There  are  a  number  of  possible  ways  to  solve  this  problem  in  3D  [36,40],  but  from  this 
point  onwards  we  focus  on  projection  images  and  a  3D  variant  thereof. 

For  two  dimensional  projection  images  in  the  x-y  planes  along  the  depth  direc¬ 
tion  (z),  we  assume  that  the  inhomogeneities  are  isolated  and  thin,  and  we  drop  the 
sum.  That  is,  we  take  the  tumor  function  to  be  zero  everywhere  except  at  the  phantom 
(object)  depth.  Thus  a  measurement,  Fourier  transform,  and  a  simple  division  in  the 
k-space  yields  the  tumor  function,  i.e 


T(p,q,zobj) 


§sc{y,q,Zd) 

AzGo{p,q,zd,z0bj) 


(9) 


where  the  subscript  “obj”  indicates  the  phantom  coordinates. 

The  inverse  Fourier  transform  of  this  quantity  enables  us  to  solve  for  absorption 
and  scattering  variations.  Using 


n  p  p-i2n (px’+qy')  f  r 

T(x',y',zobj)  =  /  /  dpdq-^j- - -  /  /  dxdye^  ^sc(x,y,zd)  ,  (10) 

JJ  AzG(p,q,zobj,z)  J  J 


we  can  solve  eq.  (2)  and  (3)  to  get  Sfia{ r)  and  Spfa (r).  Note  that  we  need  to  know  the 
background  field  at  the  object  depth  to  accurately  to  obtain  optical  properties.  Fig.(2)is 
a  flowchart  illustrating  the  algorithm  described  above. 
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Figure  2.  Simplified  flow  diagram  of  the  image  reconstruction  algorithm.  Dotted 
lines  are  used  for  optional  steps.  Brown  and  blue  indicate  real-m  filtering,  green 
indicates  G-filtering  steps. 

The  relevant  clinical  situation  is  shown  in  fig. (1) .  The  breast  is  embedded  in  a 
slab-like  box  filled  with  matching  material  such  as  Intralipid.  The  optical  properties  of 
Intralipid  are  chosen  to  be  close  to  the  average  optical  properties  of  the  breast.  Illumi¬ 
nation  comes  from  one  plane,  and  the  detectors  are  scanned  in  a  grid  on  the  opposite 
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plane.  The  amplitude  and  phase  of  the  DPSW’s  with  and  without  the  breast  are  meas¬ 
ured  and  used  to  obtain  the  scattered  DPDW.  In  our  experiments  [21,33,41]  we  also  use 
multiple  optical  wavelengths  to  obtain  spectroscopic  informational].  The  background 
field  at  different  depths  is  not  measured  and  must  be  estimated  using  numerical  models 
or  a  backpropagation  technique  [25,51].  This  is  the  main  difficulty  for  obtaining  3D 
information  in  clinical  situations  and  we  are  investigating  various  methods  to  overcome 
this  problem  [33,34]. 

2.3  Noise  amplification  at  high  spatial  frequencies 

In  eq.(9)  we  divide  two  waves  with  complex  wavenumbers  (m  and  fco).  The  scattered  field 
generally  suffers  from  noise;  the  origin  of  this  noise  can  be  experimental,  numerical,  or 
due  to  finite  sampling.  The  noise  in  the  high  spatial  frequencies  is  preferentially  amplified 
by  the  deconvolution  because  the  Green’s  function  in  the  denominator  of  equation  (9) 
decays  to  zero  faster  for  higher  values  of  p  and  q.  The  amplification  of  the  high  frequency 
noise  and  the  associated  instabilities  are  common  problems  associated  with  Fourier 
convolution  approaches  to  image  reconstruction  [25-28,44-46].  Fig. (3)  demonstrates  this 
problem. 


Figure  3.  (a)Amplitude  of  the  scattered  wave  as  a  function  of  p  for  a  fixed  q 
and  (b)  amplitude  of  the  tumor  function  plotted  in  the  k-space  for  noiseless  (only 
numerical  noise)  and  (c)  noise- added  data  from  the  single  object  (sec.  5.1)  .  The 
maximum  frequency,  |p|  =  7r/Aa:.  The  rising  “wings”  on  the  sides  are  due  to  noise. 

The  noise  effects  are  amplified  by  w  103  relative  to  the  noiseless  case  at  large  p. 

One  [37,38,41]  way  of  dealing  with  this  problem  is  to  apply  a  low-pass  spatial 
filter  in  order  to  suppress  the  high  frequency  components  of  the  tumor  function.  The 
quantitative  effects  of  this  cut  are  not  known  however  [41,45,46].  Furthermore,  the  spa¬ 
tial  resolution  of  the  images  depend  on  the  highest  k-space  frequencies  available  and 
therefore  filtering  modifies  image  resolution  (see  eq.(12))[31,41,43,44]. 

3.  Filtering  and  normalization  for  improving  depth  information  /  image 
quality 

We  have  investigated  filtering  and  normalization  techniques  which  provide  different 
types  of  information  and  lead  to  the  improvement  of  image  quality  within  the  context 
of  projection  imaging.  In  particular  we  have  found  a  robust  method  which  enables 
the  experimenter  to  estimate  the  approximate  depth  of  the  heterogeneities  below  the 
detection  plane;  the  method  does  not  require  more  data,  yet  offers  a  prescription  for 
extending  the  2D  projection  approach  to  quasi-3D  imaging  and  localization.  In  this 
section  we  briefly  outline  three  techniques  to  optimize  image  reconstruction.  In  the 
remainder  of  the  paper  we  demonstrate  their  utility.  For  simplicity  we  focus  entirely  on 
absorption  heterogeneities. 


#9187 -$15.00  US 
(C)  1999  OSA 


Received  March  01,  1999;  Revised  March  29,  1999 
12  April  1999  /  Vol.  4,  No.  8  /  OPTICS  EXPRESS  254 


3.1  Phase  only  Green’s  function  (real-m  filter) 

We  have  found  a  very  crude  filtering  technique  that  gives  fairly  accurate  positional 
information  in  two-  and  even  three-dimensional  images.  In  essence  the  method  is  based 
on  the  hypothesis  that  most  of  the  positional  information  is  contained  in  the  phase  of  our 
signals  in  k-space.  In  this  procedure  we  modify  the  Green’s  function ,G0(p,q,  zd,  z')]  in 
particular  we  put  G0(p,q,  zd,  z')  »  !*«*-*' I .  We  set  its  amplitude  to  unity,  so  only 

phase  information  in  the  Green’s  function  is  used  for  the  reconstruction.  We  also  apply 
a  Blackman  Window  on  T  (see  [43]  for  its  properties)  which  is  a  standard  windowing 
function  used  in  the  Fourier  domain  to  further  stabilise  the  image.  This  last  step  is  used 
for  better  image  quality  and  is  not  an  essential  part  of  the  algorithm.  Calculation  of 
the  projection  image  proceeds  along  similar  lines  except  for  these  steps  (see  fig.  (2)). 
Although  quantitative  information  about  optical  properties  is  largely  lost,  the  real-m 
filter  has  no  free  parameters  (except  from  the  optional  windowing  function)  and  provides 
superior  information  about  the  xyz  central  coordinates  of  heterogeneities. 

3.2  Depth  estimation  (Sj) 

We  have  empirically  found  [52]  that  the  following  set  of  operations  applied  to  the  pro¬ 
jection  images  obtained  with  the  real-m  filter  enable  us  to  accurately  estimate  the  het¬ 
erogeneity  depth  below  the  detection  surface.  Briefly,  we  first  obtain  the  tumor  function 
for  each  slice  (centered  on  some  Zj)  using  the  real-m  filter.  We  then  derive  an  image  of 
5p,a{x ,  y,  Zj)  for  each  slice,  determine  the  center  position  of  the  object  we  are  trying  to 
characterize,  and  record  the  absorption  variation  at  the  object  central  position.  Finally 
we  compute  the  quantity 


Sj 


|<5/xa  -  6fj.a\ 


axy  Spay 


(ii) 


where  the  sum  in  the  denominator  is  over  all  pixels  in  the  image  slice.  5p,a  is  the  mean 
5paxy  in  the  jth  plane.  Sj  is  in  essence  a  measure  of  the  contrast  of  the  object  in  the 
jth  plane.  We  empirically  found  that  the  value  of  Zj  for  which  Sj  is  a  maximum  closely 
approximates  the  actual  depth  of  the  object  (z0y).  This  procedure  enables  us  to  select 
a  slice  for  the  projection  image  and,  thus  provides  a  scheme  for  extending  images  to 
three  dimensions  (see  section  5.1  .1). 


3.3  High  frequency  filtering  based  on  depth  and  experimental  noise  floor  (G-filtering) 


After  we  have  obtained  this  rough  picture  of  the  object  in  the  turbid  medium,  we  repeat 
the  entire  procedure  using  the  true  Green’s  function.  From  our  discussion  in  section  2.3 
however,  it  is  apparent  that  we  need  to  devise  a  scheme  to  select  the  spatial  frequencies 
for  reconstruction.  We  carry  out  this  procedure  in  a  straightforward  way.  Note  that 
it  is  critical  for  this  procedure  that  the  scattered  waves  are  normalized  by  the  source 
strength.  For  our  convenience  we  used  a  source  strength  of  one  photon  per  unit  volume 
per  unit  time.  We  first  plot  the  raw  data  to  determine  the  experimental  noise  level.  For 
example  by  plotting  <lsc  vs.  p  and  q  (e.g.  fig. (3))  we  will  find  at  sufficiently  high  p  and 
q ,  that  4>sc  hits  a  “white”  noise  floor;  the  ratio  of  the  signal  amplitude  at  p  —  q  =  0  to 
the  signal  amplitude  at  these  high  frequencies  provides  a  measure  of  the  experimental 
signal-to-noise.  Next  we  insert  the  estimated  object  depth  z'  =  zobj  into  the  expression 
for  Go(p,  q ,  zd,  z')  ,  and  we  determine  at  what  value  of  m,  Go(p>  <?,  zd,  zobj)  drops  below  a 
threshold  based  on  our  experimental  S/N.  In  practice  we  set  this  threshold  to  be  about 
five  times  the  white  noise  floor.  We  set  f  =  0  for  all  m  greater  than  this  threshold  (see 
fig  (2)  for  the  occurance  of  this  step  in  the  general  algorithm.).  The  result  is  a  depth 
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dependent  pillbox  filter  applied  to  the  tumor  function.  The  rest  of  the  reconstruction 
is  the  same.  The  combination  technique  offers  the  possibility  of  improvement  in  the 
accuracy  of  optical  properties  and  better  images. 

At  this  point  it  may  also  be  desirable  to  recompute  Sj  using  reconstructions 
based  on  the  G-filter  technique.  Usually  the  maximum  Sj  determined  with  the  G-filter 
is  within  ^1  cm  of  the  maximum  Sj  determined  using  the  real-m  filter  (see  section  3.1 
).  In  applying  this  latter  procedure  the  high  frequency  cut-off  depends  on  depth  (zj). 

4.  Data  generation 

Our  “data”  was  derived  from  numerical  simulations.  In  particular  the  finite  difference 
method  with  partial  current  boundary  conditions  was  used  to  solve  the  diffusion  equa¬ 
tion  for  arbitrary  interior  inhomogeneities.  The  program  is  capable  of  simulating  the 
forward  problem  in  a  box  geometry.  We  used  a  point  source  with  a  strength  of  one 
photon  per  unit  volume  per  unit  time;  this  choice  of  source  term  eliminated  the  need  to 
normalize  the  scattered  wave  when  computing  the  noise  floor  and  comparing  the  thresh¬ 
old  to  the  Green’s  function.  The  source  was  modulated  at  140Mhz.  For  the  purposes  of 
this  paper  we  simulate  data  in  a  20x20x20cm  box  and  use  the  central  9x9x5cm  region 
for  “experiments” .  The  region  is  far  from  all  boundaries  and  is  a  very  good  approxima¬ 
tion  to  the  infinite  medium.  The  region  created  from  the  simulation  is  a  65  by  65  by  36 
grid  in  x,  y  and  z  with  the  coordinates  respectively,  running  from  -4.5  to  4.5cm  in  x  & 
y  and  0  to  5  cm  in  z  directions.  The  step  size  is  0.14cm  corresponding  to  a  voxel  size 
of  (0.14cm)3.  The  geometry  chosen  closely  mimics  our  clinical  set-up  [34,41].  We  use 
a  high  number  of  pixels  in  the  numerical  simulation  to  avoid  truncation  errors.  Each 
point  of  the  detection  grid  is  centered  in  the  squares  of  the  same  65x65  region  in  the 
detector  plane  at  Zj  —  Zd  =  5cm. 

Optical  properties  can  be  assigned  to  each  pixel  independently  so  that  we  are 
able  to  simulate  heterogeneities  of  arbitrary  shape.  We  use  thin  rectangular  slices  (i.e 
thin  in  z-direction)  for  objects.  Thin  slices  insure  that  our  projection  approximation  is 
reasonable.  Figs.  (4)  and  (8)  provide  3D  renderings  of  the  input  phantoms. 

Noise  is  added  to  the  data  using  a  random  number  generator  with  a  normal 
distribution  of  variance  one  and  mean  zero.  The  approach  is  tuned  to  provide  gaussian 
noise  with  variance  0.5%  for  the  amplitude  and  variance  of  0.05°  for  the  phase.  These 
probably  overestimate  our  experimental  noise.  Noise  is  added  to  heterogeneous  and 
homogeneous  measurements.  Figs.  (4)  and  (12)  show  the  amplitude  of  the  noiseless  and 
noisy  scattered  fields  obtained  in  Born  approximation. 

5.  Results 

5.1  Tissue  phantoms  with  single  slice  heterogeneity 

The  simulated  heterogeneity  is  shown  in  fig.(4).  It  has  dimensions  1.4  x  1.4  x  0.28  cm  in 
x,y  and  z  directions  respectively.  The  optical  properties  are  (ia  —  0.1  cm-1  and  \j!s  =  8 
cm-1  .  The  background  optical  properties  are  pao  =  0.02  cm-1  and  pIso  =  8  cm-1. 
The  phantom  is  centered  at  a  depth  of  2.65cm  from  the  source  plane.  It’s  transverse 
location  centers  on  y=x=  1.97cm.  The  amplitude  of  the  “measured”  scattered  field  at 
the  detector  plane  is  shown  in  its  noiseless  and  noisy  versions  in  fig. (4)  for  comparison. 
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Figure  4.  Single  Slice  Phantom:  Left  figure  shows  a  3D  rendering  of  the  phantom. 
Gray  region  has  background  properties.  The  detector  plane  is  assumed  to  be  at 
z=5cm  and  the  source  is  at  the  origin  in  z— Ocm  plane.  Amplitude  of  the  scattered 
field  in  the  detector  plane  for  the  phantom  shown  in  the  middle  (noiseless)  and  right 
(noise  added)  figures. 


5.1  .1  Real-m  imaging 

We  now  employ  the  real-m  filter  scheme  described  in  section  3.1  on  the  single  object 
data  shown  above.  Figs.  (5a)  and  (5b)  show  two  dimensional  projection  images  at  depths 
z=2.42cm  and  z=3.24cm  respectively.  We  see  that  the  object  is  clearly  apparent.  In  fig. 
(5c),  we  plot  Sj  vs  Zj\  we  see  that  the  contrast  parameter  is  peaked  at  Zj  =  2.71cm, 
and  we  use  this  value  to  make  a  “true”  real-m  filter  reconstruction  of  the  data  which  is 
shown  in  fig.(5d).  This  result  is  representative  of  all  of  the  single-object  phantoms  that 
we  tested. 


(a)  (b)  (d) 

z=  2.42  z=  3.28  z=2.71 


(c) 

Figure  5.  (a)  and  (b)  projections  at  z=2.42,  z=3.28  respectively,  (c)  Sj  vs  Zj 

(cm)  through  the  transverse  center,  peaks  at  z=2.71,  (d)  projection  at  z=2.71.  All 
with  real-m  filter. 

Although  most  of  the  information  about  optical  properties  is  lost,  the  3D  local- 
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ization  of  single  objects  is  quite  good.  The  slices  with  the  greatest  noise  are  either  near 
the  detection  or  source  planes.  The  issue  of  multiple  objects  at  different  depths  will  be 
addressed  later  in  this  paper. 

5.1  .2  G- filter  imaging 

Next  we  use  the  combination  scheme  described  in  section  3.3  to  reconstruct  the  single 
object  data.  Our  experimental  noise  threshold  was  3xl0~3  for  all  cases.  For  Zj  —  2.71cm, 
the  G-filter  was  set  for  p=q=  1.01  /cm  at  maximum,  and  the  resulting  image  in  this  plane 
is  shown  in  fig.  (6a)  .  Notice  that  the  image  is  a  little  cleaner  than  the  real-m  image  with 
most  of  the  artifacts  from  the  characteristic  ringing  near  the  corners  of  the  reconstructed 
volume.  The  optical  properties  (although  still  not  accurate  for  such  a  deep  object)  are 
more  reasonable. 
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(a)  (c) 

z=  2.7 1  z=  4.0 


Figure  6.  (a)  projection  at  z=2.71  at  real-m  peak,  (b)  Sj  vs  Zj  (cm)  through 

the  transverse  center  peaks  at  z=4cm  ,  (c)  projection  at  z=4.0.  All  obtained  with 
G-giltering. 

We  next  computed  Sj  vs  Zj  using  the  G-filter  (see  fig.  (6b)).  Notice  that  the 
peak  value  of  Sj  has  shifted  from  its  true  value.  The  G-filter  image  (at  p=q=2.36/cm) 
based  on  this  new  Zj  is  shown  in  fig.  (6c);  we  see  that  the  image  is  sharper  and  has 
better  optical  properties.  This  is  very  important  for  multiple  wavelength  images  as  we 
observe  that  the  ratio  of  reconstructed  optical  properties  is  fairly  well  approximated 
by  this  technique.  However,  its  actual  location  is  systematically  shifted  away  from  its 
true  value.  At  present  we  do  not  fully  understand  the  origin  of  this  systematic  error  in 
a  deep  way.  We  find  that  this  shift  is  related  to  the  amount  of  power  cut  out  by  the 
filter.  By  changing  the  filter  spatial  frequency  cut-off  to  a  lower  frequency  (increasing 
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the  threshold)  we  find  that  the  shift  is  increased  but  the  optical  properties  at  the  shifted 
location  however  approach  their  correct  values. 

Finally  we  note  that  since  the  maximum  spatial  frequencies  used  in  the  G-filter 
approach  depend  on  object  depth,  we  can  expect  our  resolution  to  decrease  with  increas¬ 
ing  object-detector  separation.  Following  Pattanayak  [31]  our  approximate  experimental 
resolution,  L,  depends  on  the  maximum  allowable  p  and  q,  i.e 

Pmax  +  Qmax  =  *  (12) 

In  fig.  (7)  we  show  the  change  in  resolution  with  depth  assuming  our  experimen¬ 
tal  S/N  threshold.  While  the  exact  numerical  values  depend  on  approximations  in  [31], 
the  important  qualitative  effect  to  note  is  that  the  resolution  improves  dramatically  for 
objects  near  the  detection  plane. 

15  | - 1 - 1 - 1 - 1 - 1  f  I 


1.5  2  2.5  3  3.5  4  4.5  5 

Z{cm) 

Figure  7.  Estimate  of  resolution(cm)  vs  distance  from  source  plane  (cm)  .  The 
changing  depth  dependent  cut-off  frequency  results  in  the  increase  in  resolution  with 
distance  from  source  plane  (i.e  decrease  in  resolution  with  depth  from  the  detector 
plane). 

5.2  Tissue  phantoms  with  two  slice  heterogeneities  /  three  dimensional  renderings 

We  now  simulate  two  objects  with  same  dimensions  and  optical  properties  but  in  a  more 
complicated  geometry  where  in  principle  the  scattered  waves  from  one  heterogeneity 
could  effect  the  other.  One  object  is  centered  at  3.07cm  deep  with  its  transverse  center 
at  x=y=1.97cm  and  the  other  object  is  4.35cm  deep  with  its  transverse  center  at 
x— y=-1.97cm  as  shown  in  fig. (8).  The  amplitude  of  the  “measured”  scattered  field  at 
the  detector  plane  is  shown  in  its  noiseless  and  noisy  version  in  fig.  (8)  for  comparison. 


Figure  8.  Two  slice  Phantoms:  Two  leftmost  figures  show  3D  renderings.  Gray 
region  has  background  properties.  The  detector  plane  is  at  z=5cm  and  the  source  is 
at  the  origin  in  z=0cm  plane.  Amplitude  of  the  scattered  field  at  the  detector  plane 
is  shown  in  the  two  rightmost  figures.  The  left  is  the  noiseless  data  and  the  right 
shows  the  data  after  adding  noise. 

We  performed  essentially  the  same  set  of  reconstruction  procedures  as  described 
in  sections  3.  and  5.1  .  In  fig. (9a)  we  show  that  Sj  has  two  maxima  in  different  planes 
depending  on  the  heterogeneity  under  consideration  using  real-m  filter.  Note  that  the 
real-m  images  were  more  noisy  by  comparison  to  the  single  object  phantom.  We  again 
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find  that  the  optical  properties  are  most  accurate  where  Sj  has  its  shifted  maxima  with 
the  G-filter.  These  images  are  shown  in  Fig  (9b, c).  Absolute  values  are  closer  to  the  real 
values  for  shallow  objects  and  the  positional  shift  also  depends  on  the  depth. 


y(cm) 


.2  0  2  4  -4  -2  0  2  4  x(cm) 

x(cm) 

(b)  (c) 

z=  3.71  z=  4.42 

Figure  9.  (a)  Circles  (crosses)  show  Sj  vs  Zj  (cm)  through  the  center  region  of 

deeper  (shallower)  object  obtained  with  m-filter.  Peak  for  both  objects  are  exact. 

Then  with  G-filter  we  get  projections  at  (b)  z=3.71  and  (c)  at  z=4.41. 


Finally,  by  calculating  Sj  at  all  the  voxels  we  can  generate  three  dimensional 
images  of  the  entire  reconstructed  volume.  That  is,  for  each  of  the  36  projections  we 
calculate  the  value  of  Sj  for  all  transverse  pixels  (65x65=  4225  values  in  each  projection) 
and  plot  all  S/s  in  three  dimensions.  In  fig. (10a)  and  (10b)  we  show  an  isosurface 
rendering  at  value  Sj  =  0.042  at  two  viewing  angles.  Two  dimensional  projections  in 
plane  y=2cm  and  plane  y=-2cm  are  also  shown  in  fig.  (10c)  and  fig.(10d)  to  provide  a 
better  sense  of  contrast.  Comparing  this  to  fig. (8)  we  see  that  the  shallow  object  is  very 
well  reconstructed  both  in  terms  of  shape  and  location.  The  deeper  object  has  lower 
resolution  and  is  shifted  from  its  original  location.  In  the  corners  we  see  the  characteristic 
ringing  effects  from  Fourier  domain  reconstructions.  These  latter  artifacts,  however,  are 
usually  fairly  easy  to  identify.  We  find  that,  even  when  the  objects  are  close  to  the  image 
edges  we  are  capable  of  separating  objects  from  ringing  artifacts.  We  are  investigating 
presently  whether  it  is  possible  to  extract  any  quantitative  or  qualitative  optical  property 
information  from  these  images. 
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Figure  10.  Three  Dimensional  Rendering  of  the  G-filter  reconstruction.  An  iso¬ 
surface  at  Sj  =  0.042  is  shown  in  two  different  angles  in  (a)  and  (b).  In  (c)  and 
(d)  images  of  Sj  in  x-z  plane  through  y=2cm  and  y=-2cm  are  shown  respectively. 
Compare  the  results  to  that  of  fig.  (8) 
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6.  Conclusion 


We  have  studied  projection  images  of  turbid  media  based  on  the  angular-spectrum 
representation.  We  demonstrated  two  filtering  schemes.  One  of  them,  the  real-m  fil¬ 
ter,  provided  object  location  information  with  no  free  parameters;  the  other  method, 
the  G-filter,  introduced  a  systematic  approach  for  eliminating  high  spatial  frequencies 
from  the  reconstructions.  Image  localization  in  3D  was  demonstrated  using  both  filters 
by  maximizing  the  object  signal-to-noise  ( Sj )  on  a  slice-by-slice  basis.  Despite  these 
successes,  some  image  artifacts  remain  and  must  be  addressed  in  future  work.  These 
include:  noise  near  the  detection  plane  and  image  boundaries,  systematic  shifts  of  the 
object  location  depending  on  filtering  schemes  and  noise  floors,  and  optical  property 
accuracy. 
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Abstract 


An  angular  spectrum  algorithm  is  presented  for  fast,  near-field  diffraction 
tomographic  ilmaging  with  diffuse  photon  density  waves  in  highly  scattering 
media.  A!  general  relation  in  K-space  is  derived  that  connects  the  spatial  vari¬ 
ations  of  the  optical  properties  of  heterogeneities  to  the  spatial  spectra  of  the 
measured  scattered  diffuse  photon  density  waves.  The  theory  is  verified  exper¬ 
imentally  for  situations  when  boundary  effects  can  be  neglected.  We  further 
describe  Kowi  to  incorporate  boundary  conditions  into  this  angular  spectrum 
algorithm  for  a  turbid  medium  of  finite  size,  the  slab  medium.  Limitations 
and  potential  improvements  of  the  near-field  diffraction  tomography  are  also 
discussed. 
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I.  INTRODUCTION 


Optical  radiation  was  used  to  image  breast  tumors  by  the  shadowing  effect  as  early  as  the 
1920’s  [1],  However,  recent  advances  in  light  generation  andl  detection,  along  with  improve¬ 
ments  in  our  theoretical  understanding  of  near,  infrared  flNIR)  light  propagation  in  tissue-like 
highly  scattering  turbid'  media  have  opened  new  possibilities  for  optical  imaging  of;  the  inte¬ 
rior  of  thick  biological  tissues  [2].  In  highly  scattering  media  such  as  biological  tissue,  light 
propagation  is  described  adequately  within  the  diffusion  model  of  photon  transport  [3-5] . 
It  has  been  shown  by  several  investigators  that!  diffuse  ph'oton  density  waves,  which'  are  cre¬ 
ated  inside  highly;  scattering  media  by  an  intensity  modulated  light  source,  obey  a  Helmholtz 
wave  equation  with  a  complex  wave  number  [6,7].  In  spite  of  the  complexities  resulting  from 
strong  tissue  scattering,  the  new  modality  with  diffusing  photons  offers  many  attractive  fea¬ 
tures  for  imaging  thick  tissue.  These  features  include  non-invasiveness,  low  cost,  and  unique 
optical  contrast  and  spectroscopic  signatures  with  clinical  and  physiological  relevance  [8,9]. 

The  goal  of  optical  imaging  with  diffusive  photon  density  waves  is  to  reconstruct  a  low 
resolution  map  of  the  heterogeneous  absorption  and  scattering  properties  from  the  measure¬ 
ments  of  diffuse  photons  on  the  sample  surface  or  within  the  sample.  Image  reconstruction 
entails  solving  the  inverse  problem.  Most  quantitative  optical  image  reconstruction  algo¬ 
rithms  such  as  the  Algebraic  Reconstruction  Technique  (ART)  and  the  Spontaneous  Iterative 
Reconstruction  Technique  (SIRT)  [110],  the  Newton-Riaphson  technique  combined  with  finite 
element  numerical  method  [111-13],  the  Conjugate  Gradient  Descent  technique  [14],  and  the 
Singular  Value  Decomposition  (SVD)  [15],  rely  on  iterative  schemes  in  a  least-square  sense. 
The  optical  image  reconstruction  therefore  requires  a  significant  amount  of  computational 
resources  and  time. 

Recently;  we  showed  that  by  essentially  following  the  techniques  of  diffraction  tomog¬ 
raphy  |16,17],  it  is  possible  to  rapidly  reconstruct  thin  slice  and  spherical  objects  whose 
absorption  and/or  scattering  parameters  differ  from  the  background  homogeneous  scatter¬ 
ing  medium  |18].  The  image  reconstruction  algorithm  employing  the  diffraction  tomography 
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technique  (called  angular  spectrum  algorithm!  ini  this  paper]  is  rapid,  permitting  object  local¬ 
ization  and  characterization  in  ~  1000  vodume-elemiemt  samples  on  sub-second  computational 
time  scales.  The  primary  purpose  of  this  paper  is  to  provide  a  more  complete  discussion  of: 
the  results  reported  in  those  earlier  papers.  In  addition,  we  will  provide  detailed  analysis  of 
this  algorithm  Incorporating  the  effects  of  finite  boundaries.  We  will  first  derive  the  general 
integral  solution  of  the  total  and  scattered  photon  density  waves  in  a  heterogeneous  turbid 
medium  within  the  first  order  Born  approximation  (Section  II,  III  and  IV).  We  next  derive 
a  relation  in  Iv-space  between  the  spatial  spectrum  of  the  heterogeneity  function  and  th'e 
spatial  spectrum  of  the  measured  scattered  diffuse  photon  density  wave  (Section  V  (A))  .  Ex¬ 
perimental  results  are  presented  to  verify  the  feasibility  of  the  angular  spectrum  algorithm 
for  image  reconstruction.  Some  limitations  and  potential  improvements  of  the  diffraction 
tomography  are  also  discussed  (Section  VI).  Finally  we  present  a  discussion  about  how  to 
incorporate  boundary  conditions  into  this  angular  spectrum  algorithm  for  a  turbid  medium 
of  finite  size,  in  particular,  the  slab  medium  and  the  semi-infinite  medium  (Section  VII). 

II.  PHOTON  DIFFUSION  EQUATION  IN  HETEROGENEOUS  MEDIA  -  A 

PERTURBATION  APPROACH 

Light  transport,  in  highly  scattering  turbid  media  is  often  well  described  by  photon  diffu¬ 
sion  |2].  Consider  a  light,  source  at  r,  with  its  intensity  sinusoidally  modulated  at  modulation 
frequency  /.  e.g.,  the  source  term  is  S(r,  t)  —  S{v)e~luJ  1  =  Moe~‘UjJt5(r  —  rs),  where  lo  =  2tt  f 
is  the  angular  source  modulation  frequency1,  M0  is  the  source  strength  representing  the 
number  of  photons  emitted  per  second.  Consider  steady  state  photon  diffusion  in  which  the 
photon  fluemce  <fi(]n.  £]  has  the  same  time  dependence  as  the  source,  i.e.,  $(r,£)  -  $(r)e-*w*. 
Tt.  is  straight,  forward  to  show  that  the  photon  fluence  <&(r)  satisfies  the  photon  diffusion 


1  The  continuous- wave  (OW)  case  is  a  special  case  where  u>  —  0  and  the  frequency  domain  analysis 


nail!  be  readily  applied  to  th'e  OW  case. 
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equation  [3-3]: 


V  •  (D  V<$(r])  -  u^aH r)  +'  =  -vS{ r)  .  (1) 

Here  the  common  time  dependence  exp(-i  u)  t)  ofl  the  fluenee  <B(r)  and  the  source  .9 (in)  are 
omitted,  v  is  the  speed  of  light  in  the  turbid  medium;  D  =  is  photon  diffusion:  coefficient j 
fia  and  /x*  are  respectively  the  optical  absorption  and  reduced  scattering  coefficients.  Since 
V  •  (DW) |  =  +  SAD  -  V$  and  S7D  =  =  -^r  ^rS  Eq.  (]1)  can  be  rewritten  as 

V24(r)  41 3^(~i,,**+“J)$(r.)  -  M  .  V4(r)  =  -3,4  S(r)  .  (2) 

V  Ms 

Tm  a  homogeneous  medium,  the  absorption  and  scattering  coefficients  (jiM  and  p'.0)  are 
constant,  and  this  above  equation  reduces  to  a  simple  Helmholtz  equation: 

(V2  +  *„2)<Mr)  =  —3,4  Sir)  .  (3) 

Here  the  wa\ie  number  D  is  complex  and  Uq  =  [3 /jfs ( - )] 1//2  with  Irn(ko)  >  0  to  ensure 
that  the  photon  density  goes  to  zero  at  a  large  distance.  The  simplest  solution  of  this 
Helmholtz  equation  shows  that  photons  propagate  collectively  as  a  damped  spherical  wave 
outward  from  the  source,  i.e.,  <ho(r)  oc  eX^|t-rJ|~- 

Tn  am  optically  heterogeneous  turbid  medium,  the  spherical  wave-fronts  of  the  homoge¬ 
neous  incident  wave  are  distorted  by  inhomogeneities  and  scattered  waves  are  produced. 
The  concept  of  the  scattered  wave  is  illustrated  in  Fig.  1.  The  total  photon  density  wave 
<R(]r)|  is  the  sum  of  the  incident  wave  <T>0(r)  and  the  scattered  wave  $sc(r): 

$(r.)  =  <Il0(n)  +  $sc(r)  .  (4) 

This  formalism  is  not  limited  to  an  infinite  geometry.  In  general  the  incident  wave  $0(>*) 
corresponds  to  the  photon  density  waive  in  a  homogeneous  turbid  medium  for  an  arbitrary 
geometry,  and  wie  interchangeably  use  another  name  -  "background  photon  density  wave”  - 
for  4>0(iri)i  the  scattered  wave  4>„.(:n)  represents  the  perturbation  of'  the  incident  wave  in  the 
presence  of  imhnmogenieities  for  the  same  geometry.  The  scattered  wave  is  determined  by 
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characteristics  of  thle  inhomogeneity  such'  as  its  size,  shape,  position  and  its  absorption  and 
scattering  properties.  The  goal  of  optical  tomography  is  to  reconstruct  the  map  of  these 
heterogeneous  optical  properties  from  measurements  of  the  scattered  wave.  We  first  derive 
thle  photon  diffusion  equation  obeyed  by  the  total  photon  density  wave  within  the  first  order 
Bonn  Approximation. 

In  a  heterogeneous  medium  we  write  the  optical  properties  (/-«„  and:  p',)  as  the  sum  of 
background  optical  properties  (^ao;  n'so)  an(l  the  variations  relative  to  the  background 

W.  ie- 

Ha  =  Hoi)  +  *Ha  ,  (|S) 

h's  =  h'so  +  Sh's  ■  (6) 


We  consider  the  case  of  weak  optical  inhomogeneities  where  S/ia  <C  Ha o  and  5p!s  <C  /z'0.  The 
optical  inhomogeneilties  introduce  a  weak  perturbation  to  the  incident  wave,  i.e.,  |$sc(r)|  <C 
|<f»0(n)|.  Substituting  Eqs.  (5)  and  (6)  into  Eq.  (2)  and  keeping  only  the  zeroth  and  first 
order  terms  in  optical  property  variations  (i.e.  S/ia,  5/j!s)  as  well  as  in  the  scattered  wave 
wie  find 

(V5  +,  ti)i »(r)  =  ~fl  +  yp-]  S( r)  -  TU r)  -  T„( r)  ,  (7) 

uo  Hso 

where  we  have  introduced  the  heterogeneity  functions  Tabs(.r)  and  Tsc(r)  representing  the 
perturbations  due  to  the  absorption  and  scattering  variations.  They  are 


Tahsi,r)  =  -T:0(r)  Stia(r)  , 

uo 


_  .  ,  3D0/:f( 


*oW  V.W  -  ■  V<I>o(r)  . 


n'so 


(8) 

(9) 

Here  D0  =  ^4-  is  the  background  photon  diffusion  coefficient;  /c0  =  [{—v/.ia 0  +  i  lS)/D$'I2 
is  the  wavenumber  of  the  incident  diffuse  photon  density  wave  <f>0(r).  Note  that  ^4^-  S( r) 
in  the  first  term  on  the  right-hand  side  of  Bq.  (7)  is  zero  as  long  as  the  source  is  outside 
the  inhomogeneity  (which  is  generally  the  case  in  practice),  and  therefore  we  can  drop  this 
term  from  Fkj.  (7).  Tm  addition  we  assume,  for  simplicity,  that:  the  scattering  varies  slowly  in 
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space  and  the  second  terra  In  Eq.  (9)  vMo+'Wr)1  •  V$0(in)  can  then  he  neglected.  We  thus 
have  the  following  simplified  equation  within  the  first  order  Born  Approximation) 

(V?  +  kl)mr)  =  ~  S1(r)  -  T.( r)  ,  (10) 

where  T(r)  =  Ta6s(jr)  +  Tsc(r)  and  the  heterogeneity  functions  Toi)S(jr)  and  Tsc(r)  are  given 
by 

litoW  = -s-®oW  i«.w  (11) 

M) 

T>e(r)  =  !M1,0(r)V,(r)  (112) 

We  see  that  the  heterogeneity  funotions  due  to  the  optical  inhomogeneities  can  be  treated 
as  equivalent,  "source”  terms  which  give  rise  to  the  scattered  component  <3>.5C(r)  of  the  total 
diffuse  photon  density  wave  T(r). 


III.  TOTAL  DIFFUSE  PHOTON  DENSITY  WAVE  IN  HETEROGENEOUS 
TURBID  MEDIA  —  THE  GREEN’S  FUNCTION  APPROACH 


We  will  take  a  Green’s  function  approach  to  derive  the  total  and  therefore  the  scattered 
photon  density  wave  in  a.  heterogeneous  highly  scattering  medium.  Consider  the  Green’s 
function  in  turbid  media  which  satisfies 


(V2  +  q)G(|r,r‘)  =  -i(r,r‘),  (13) 

where  k0  —  [(|— w/irt0  4-  i  oj)/  D0]l/2  i:s  the  wavenumber  of!  the  inoident  photon  density  wave 
(the  same  as  in  Flq.  (10)')'.  Using  the  Green’s  theorem,  we  obtain  an  integral  expression  for 
the  total  optical  density  wave  <b(r): 


n]  =  ^  J1  S(r')  G(r ,  r')  dV  -H  /  T(jr')  G(r,  r')  dV 


v 


V 


(14) 


The  first,  term  on  the  right-handi  side  of  Eq.  (14)  is  a  volume  integral  of  the  light  source 
over  the  entire  turhid  medium.  Tt  gives  us  the  incident  wave.  The  second  term  is  a  volume 
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integral  of!  the  heterogeneity  function  ovier  the  entire  turbid  medium  and  it  determines  the 
perturbation  resulting  from  the  optical  heterogeneities.  The  third  term  is  a  surface  integral 
over  the  closed  surface  of  the  entire  turbid  medium.  It  takes  into  account  the  boundary 
effects  on  the  total  photon  density  wave,  and  it  includes  contributions  to  the  total  photon 
density  wave  from  both  the  incident  wave  and  the  scattered  wave  on  the  boundary,  n'  in  the 
surface  term  denotes  the  surface  normal  pointing  outward.  For  an  infinite  heterogeneous 
medium,  this  surface  term  is  zero  since  the  enclosure  surface  of  an  infinite  medium  is  at 
infinity.  Therefore  the  scattered  wave  can  be  simply  separated  from  the  incident  wave. 
For  a  finite  turbid  medium,  however  the  separation  of  the  incident  wave  component  from 
the  scattered  wave  component  in  the  surface  term  is  generally  difficult.  It  is  advantageous 
therefore  to  remove  the  surface  integral  from  the  total  photon  density  wave  by  choosing  an 
appropriate  Green’s  function.  We  will  consider  this  complicated  (yet  more  realistic)  case  at 
the  end  of  this  paper.  We  will  start  with  a  simple  case  -  the  infinite  geometry  case. 


IV.  SCATTERED  WAVE  IN  INFINITE  HETEROGENEOUS  TURBID  MEDIA 


Als  shown  in  Fig.  1 ,  in  the  presence  of  optical  heterogeneities,  the  total  photon  density 
wave  consists  of  the  incident  wave  and  the  scattered  wave,  and  the  scattered  wave  carries 
the  information  ofi  the  optical  inhomogeneities.  Starting  with  Eq.  (14),  we  can  easily  find 
the  relation  between  the  scattered  wave  $sc(r)  and  optical  inhomogeneities,  for  example 
through  the  heterogeneity  function  Tabx(r)  and  Tsc(r).  For  an  infinite  geometry,  the  surface 
integral  in  Eq,  (14)  disappears.  The  incident  wave  in  this  case  is  given  by  the  first  term 
d volume  integral  of  the  source)  on  the  niight-hand  side  of  Eq.  (.14).  For  an  infinite  geometry, 
the  Green’s  function  t7o(]n,n')  which  satisfies  Eq.  (13)  is  [19] 


exppoln  -  r'l) 


(15) 


47r|r  —  r'| 

this  Green’s  function  and  considering  an  infinite  medium  with  a  point  source  at  rs, 
i.e.,  .9(|n'|  =  A/odflr'  -  na).  we  can  readily  obtain  the  incident  wave,  i.e. 


where  Ma  is  number  of  photons  emitted  from  the  source  per  second.  We  see  that  this 
incident  wave  iis  of  course  the  same  as  the  one  obtained  by  directly  solving  the  Helmholtz 
Equation.  The  scattered  wave  (by  definition)  is  the  difference  between  the  total  photon 
density  wave  <f>(r)  and  the  incident  wave  <t>0(r),  which  is  related  to  the  volume  integral  of 
the  heterogeneity  function,  ii.e. 


Mr)  =  M)  -  <Mr)  =  J\T(r'}  <2o(r,r'J  dV 

V 


(117) 


V.  IMAGE  RECONSTRUCTION  ALGORITHM  AND  EXPERIMENTAL 

RESULTS 

The  scattered  wave  depends  on  the  heterogeneity  function.  In  practice  the  scattered 
wave  can  be  obtained  from  measurements  and  the  knowledge  of  the  incident  wave.  Given 
the  scattered  wave,  how  can  one  obtain  the  heterogeneity  function  and  thus  S/ua(r)  and 
di/i' (]n)|?.  The  approach  we  take  here  employs  the  angular  spectrum  analysis  of  the  scattered 
wave.  Tn  this  approach  we  relate  the  spatial  spectrum  of  the  scattered  wave  to  the  spatial 
spectrum  of  the  heterogeneity  function.  The  basic  ideas  are  similar  to  those  of  diffraction 
tomography  [16].  The  analysis  involves  forward  and!  inverse  Fourier  transforms  following  the 
convent  inns  given  in  Appendix  A. 


A.  The  Angular  Spectrum  Algorithm 

The  detection  geometry  we  consider  for  the  angular  spectrum  algorithm  is  a  2-D  planar 
geometry.  As  shown  in  Fig.  2.  the  scattered  wave  $sc(r)  is  determined  at  a  plane  z  =  zd 
from  a  set  of  measurements  in  that,  plane.  Eq.  (17)  tells  us  that  the  scattered  wave  T5C(r) 
is  the  convolution  of  the  heterogeneity  function  T.(r)  with  thle  Green’s  function  Gq (r,  E) .  In 
order  to  reveal  the  relation  between  the  scattered  wave  and  the  heterogeneity  function  in 
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K-space,  we  first,  expand  the  Green’s  function  in  terms  of  plane  warns  in  two-dimensions, 


G„(  r„,  r')  =  IJ  dpdq  G„(ft  9,  %,  **) 

“OC 
*f  oc 

=  JI  dpdq  ) 


where  (Jp,  (?)  are  the  2-D  spatial  frequencies  with  respect  to  the  x-y  coordinates.  In  the 
second  line  of  the  above  equation,  we  have  employed  the  Wey.1  expansion  of  the  Green’s 
function  [20],  i.e. 


where  m  =  [h?  -  (27r)2(jp2  4-  q1))^2  and  Im(m)  >  0.  The  derivation  of,  the  Weyl  expansion 
of  the  Green’s  function  is  given  in  Appendix  B. 

Note  fliat.  Eq.  (18)  is  the  angular  spectrum  representation  of  the  Green’s  function,  a 
solution  of  the  wave  equation  with  a  point  source  at  (x\y\z')  (see  Eq.  (15)).  At  any  point 
inside  the  halfi  space  to  the  right  (or  left)  of  the  source,  there  are  eigen-plane  waves  in 
fife  x-y  plane  whose  amplitudes  and  phases  vary  with  the  distance  from  the  source  \zd  -  z'\. 
Because  of;  the  large  positive  imaginary  part  of  m.  the  amplitude  decays  exponentially  versus 
the  perpendicular  distance  \zd  -  z'\  away  from  the  source  point.  Plane  waves  with  large 
spatial  frequencies  (p.  q)  (and  therefore  a  large  imaginary  part  of  m)  will  have  negligible 
amplitudes.  This  is  the  characteristic,  difference  between  diffuse  photon  density  waves  and 
ordinary  diffractive  electromagnetic  waves  in  lossless  dielectric  media.  These  plane  waves  will 
lie  scattered  by  optical  inhomogeneities  and  their  resulting  amplitudes  and  phases  will  carry 
information  about  the  absorption  and/or  scattering  characteristics  of  the  inhomogeneities. 

Tf  we  substitute  the  angular  spectrum  representation  of  the  Green’s  function  (Eq.  (18)) 
into  the  volume  integral  of  the  scattered  wave  given  by  Eq.  (17),  after  simple  algebraic 
manipulation  and  interchanging  the  order  of  integrations,  we  obtain  the  following  represen¬ 
tation,  known  as  the  angular  spmtmum,  representa  tion  of  the  scattered  wave: 


f) 


+0C 

$«:(**)  =  JJ  dpdq  e~i7x{md^qVd)  J  dlz!  G0(p,  <7,  Ab,  A  "F(p,  <?,  d]  ,  (20] 

—  OO 

where  T(p,  q,zr)  is  the  2-D  spatial  spectirumi  (Fourier  transform)!  of  the  heterogeneity  func¬ 
tion,  i.e. 

-foe 

f  (p,  q ,  z')  =  f\j\  dx'dy'  T(x',  y',  z)  ei2r  te'+w'l  .  (21) 

—  OO 

Taking  the  2-D  Fourier  transform  of  the  scattered  wave  <h4C(rj)  in  the  detection  plane  at 
21  =  i.e., 

-foe 

*«M  fj  <l»k  *«(p, ft  zd)  e-°“ in+m)  ,  (22) 

—  OO 

and  comparing  Eq.  (20)  and  Eq.  (22),  we  then  obtain  the  relation  between  the  spatial 
spectrum  of  the  scattered  wave  and  the  spatial  spectrum  of  the  heterogeneity  function  at 
any  given  spatial  frequency  (p,q),  i.e. 

OC 

$«:(?>,  q„  zd)  =  J  dz  Gq(p,  q,  zd,  z1)  T(p,  q,  z)  .  (23) 

—  OO 

Without,  losing  generality,  we  assume  the  optical  heterogeneities  are  below  the  detection 
plane.  This  assumption  enables  us  to  remove  the  absolute  value  sign  in  the  Weyl  expansion 
in  Eq.  (19]  since  a,;  —  z1  >  0.  Wb  also  assume  the  heterogeneities  are  localized  between  the 
detection  plane  at  z  =  zd  and  a  plane  at  z  —  z{}.  Thus  we  need  consider  only  the  interval 
between  (71=20,  z=zto]  for  the  integral  in  Eq.  (23).  Dividing  the  turbid  medium  between  the 
plane  at  2  =  zq  and  the  detection  plane  into  slices,  we  can  rewrite  Eq.  (23)  in  the  following 
form  of  discretized  summation  1 


dwfb,  q,  zd)  =  )T  A 2  T (p,  q,  %]  G0(p,  q,  zd,  zj ) 

.7=1 

ji=  l 


(24) 


where  in  the  second  line  we  substitute  the  Green’s  function  <70(p,  q,  zd,  Zj)  with  its  Weyl 
expansion  (Eq.  (19))  ;  Azi  is  the  discretized  step  size  along  the  z-direction  and  N-  is  the  total 
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number  of!  slices  in  the  z-dirediiom.  Ideally  the  discretization  step,  size  A  a  meeds  to  be  as 
small  as  possible.  In  practice  we  choose  A 2a  to  be  a  few  random!  walk  steps  (i.e.  ^r). 

Eq.  ()24)  implies  that  at  any  given  spatial  frequency  Op,  q'h  the  heterogeneity,  functions  at 
different  depth1  z/s  can  be  thought  of  as  the  “source  terras”  for  the  scattered  wave.  The  plane 
waves  arising  from  different  slices  propagate  along  the  z  direction  to  the  detection  plane. 
During]  tih'e  propagation  these  plane  waves  experience  different  amplitude  attenuation  and 
phase  shifts  which  are  given  by  el"^Zd~z^ /nv,  where  m  —  [Aijj  —  (27r)2('p2  +  g2)]l'/2  is  a  complex 
number  with  lm(rn)  >  0;  the  scattered  wave  detected  at  the  detection  plane  z  -  zd  is  thus  a 
sum  of:  plane  waves  originating  from  the  heterogeneity  functions  at  different  depths.  In  Fig.  3 
we  illustrate  this  concept.  I11  this  figure  we  consider  two  non-zero  heterogeneity  functions 
f|  (p,.  a,)  and  To(p.  z2)  corresponding  to  plane  waves  along  the  x-direction  in  the  x-z  plane 
fli.e..  y  =  0]  with  a  spatial  frequency  p  at  depth  z}  and  z2.  The  perturbations  from  these 
two  slices  propagate  to  the  detection  plane  with  a  phase  shift  and  amplitude  attenuation 
factor  /m.  Alt  the  detection  plane  the  perturbations  from  these  two  slices  add  up 

to  miake  a  scattered  wave  <hs0(p,  zdJ  at  the  same  spatial  frequency  p. 

Tni  K-space  tlie  propagation  of  the  perturbation  T(p,q,Zj )  at  different  depths  zd  -  Zj  is 
weighted  by  the  amplitude  attenuation  and  phase  shift  given  by  the  Weyl  expansion  of  the 
Qreem’s  function  (roflft,  q,  zd,  Zj)  =  frimN-^)/(2m).  Recall  m  =  [&g  -  (2tt)2(p2  +.?2)],/2  with 
Tm(m)  >  0.  therefore  the  amplitude  and:  phase  of  the  Weyl  expansion  Ga(p,q,  zd,  Zj)  depend 
on  the  spatial  frequency,  (ft,  <7)  at:  a  given  depth  zd  -  zj.  The  amplitude  decays  more  quickly 
as  the  spatial  frequencies  (]p,  q]  increase,  and  the  Green’s  function  effectively  acts  as  a  low 
pass  filter  in  K-space. 

For  spatial  frequencies  (ft.  q)  with  the  range  of  (0,  1.6)  cm-1,  we  plot  the  amplitude  and 
phase  of  the  Weyl  expansion  (  ~  in  Dig.  4  (a,  b)  assuming  the  depth  is  zd-Zj  = 

1  cm.  Tn  calculating  the  incident,  diffuse  wave  wavenumber  A:0  =  [(-u  pao  +  i  u)/D 0]1/2  we 
choose  background  optical  properties  pa(\  =  0.02  cm-1  and  p'sQ  =  8.0  cm-1,  and  a  140 
MHz  modulation  frequency.  The  resultant  wavenumber  is  | Ay |  ~  1.1  cm-1.  We  find  that 
the  amplitude  attenuates  by  -w  orders  of  magnitude  when  the  spatial  frequencies  (p,  q) 
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increase  from  (0,  0)  cm-1  to  (1.6,  1.6)  cm-1.  In  practice  the  maximum  spatial  frequency 
is  determined  by  the  Nyquist  sampling  frequency,  i.e.  qmfUX  =  l/2Ax  sa  0.833  cm-1  for  a 
scanning  step,  size  Ax  —  0.6  cm.  In  Flig.  4  (c,  d),  we  also  plot  the  amplitude  attenuation 
and  phase  shift  versus  depth  for  given  spatial  frequencies,  i.e.  (0.1,  0.1)  era-1  and1  (0.-3, 
0.5)  cm-1.  The  amplitude  attenuates  exponentially;  and  the  phase  shift  increases  linearly 
as  we  consider  the  perturbation  from  deeper  slices.  Again  as  already  shown  in  Fig.  4  (c.), 
the  amplitude  attenuates  much  faster  at  spatial  frequencies  (0.5,  0.9)  cm-1  than  at  (0.1, 
0.1)  cm-1.  At  any  given  depth!  (jzd  -  Zj),  those  plane  waves  with  sufficiently  large  spatial 
frequencies  (p,  q]  have  negligible  contribution  to  the  scattered  wave,  and  therefore  carry  less 
information  about  the  inhomogeneities. 


B.  2-D  Projection  Imaging 


2-D  photographic  images  have  been  used  by  radiologists  for  many  years.  In  order  to 
acquire  2-D  photographic- type  projection  images,  we  make  a  “thin'’  slice  approximation  by 
replacing  a ,•  on  the  left  hand  side  of  Eq.  (24)  with  the  estimated  slice  position  of  the  object. 
We  then  drop  the  sum  over  all  other  z/s  and  obtain  the  following  simple  relation  at  any  given 
spatial  frequency  (jp,  q)  in  K-space  between  the  heterogeneity  function  at  depth  z  =  zobj  and 
the  measured  scattered  wave  at  the  detection  plane  z  —  z&\ 


T(p.  q,  zobj) 


4c(?>,  <h  Zd) 
AzG()<lp,q,zd,zobj) 


2m, 

JMz 


<lso(p,  q,,  zd)  e 


- im(zd-zobj ) 


(25) 


This  “thin”  slice  approximation  may  be  adequate  since  we  are  often  interested  in  early 
tumors  whose  size  wiill  he  of  the  order  of  slice  thickness  of  ~0.5  cm,  and  thus  can  be 
considered  thin.  As  we  discussed  at  the  end'  of  Section  V(A),  plane  waves  in  K-space 
with  large  spatial  frequencies  (p.q)  are  attenuated  quickly  as  they  propagate  within  the 
turbid  media.  The  largest  detectable  spatial  frequencies  are  determined  by  the  sensitivity 
and  si gual-to- noise  ratio  of  the  defection  system.  Fbr  image  reconstruction,  we  neglect  the 


12 


high  spatial  frequently  components  (p.  q)  ini  tilie  heterogeneity,  function  T(p,  q.  2U,')  when 
Im(rn )  >  3.5 Jm(/c0),  by  using]  an  empirical  "nn-cut”  filter  in  K-space  [21], 

When  the  heterogeneity  function  in  K-sp;ace,  T(p,  q.  zobj).  is  determined  by  Eq.  (25), 
we  can  then  take  the  inverse  2-D  Fourier  transform  of  T(\p,  q,  zaHj)  to  obtain  the  tumor 
function  Ti(z,y,zohj)  in  the  real  x-y  space  at  the  depth  of:  the  heterogeneity  2  =  zohj. 
We  derive  a  2-D  photographic  image  of  the  optical  properties  using  Eqs.  (11)  and  (12): 
for  example,  8i.ia(x,y,zobj)  =  Tah,(x,y,zobj)/(-^  ®0(x,y,zobj))  for  absorbing  objects,  and 
6(/J:x/y,zohj)  =  Tscfo,y,zobj)/(^± l  %(x,y,zobj))  for  scattering  objects.  Note  that  for 
a  purely  absorbing  or  scattering  object,  either  a  frequency  domain  (modulation  frequency 
/  ^  0)  or  a  continuous-wave  (iCW,  /  =  0)  DPDW  can  be  employed  to  extract  the  absorp¬ 
tion  or  scattering  variation;  but  for  objects  having  both  absorption  and  scattering  variations, 
a  GW  D'PDW.  is  not  sufficient  to  separate  the  absorption  and  scattering.  Details  will  be 
discussed  in  Section  VI(A). 

Consider  next  a  oase  where  the  optical  heterogeneities  are  located  within  a  ‘‘thin”  slice 
at.  a  =  ziohj  (see  Big.  5).  If  the.  slice  thickness  Az  is  less  than  a  few  transport  mean  free  path- 
lemgtbs  fl/;(/4o  +  //.ft0j],  the  heterogeneity  function  within  this  thin  slice  is  approximately 
uniform,  therefore  Eq.  (25)  provides  a  quite  accurate  relation  between  the  heterogeneity 
function  and  tilie  scattered  wave  in  K-space,  and  optical  properties  of  the  heterogeneity  can 
further  be  deduced  quite  accurately.  For  thicker  objects  (i.e.  thickness  >4  mm),  the  average 
over  the  size  ofi  the  object  weighted  by  the  sum  of  exponential  amplitude  and  phase  factors 
jm  prOTides  only  an  approximate  relation  between  the  heterogeneity  function  and 
the  scattered  wave.  However  we;  find  that  the  relative  optical  properties  of  multiple  objects 
can  still  be  reconstructed  with  an  reasonable  accuracy. 

Obviously  the  image  reconstruction  involves  only  2-D  forward  and  inverse  Fourier  trans¬ 
forms,  and  no  iterative  schemes  are  needed;  therefore  this  angular  spectrum  algorithm  is 
very  rapid. 
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C.  a  priori  Depth'  Information  and  Perspectives  of  3-D  Imaging 


From  the  derivation  we  notice  that  in  principle,  this  K-space  spectrum  analysis  algorithm 
should  work  well  when  the  optical  heterogeneities  are  confined  wiitblim  a  thin  slice.  The 
reconstruction  then  provides  a  2-D  photographic  projection  image  of  the  optical  properties 
given  a  priori  information  about  the  depth  of  the  heterogeneity.  Since  the  heterogeneity 
function  dtherefore  the  optical  properties  of  the  heterogeneities)!  is  related  to  the  scattered 
wave  via  the  Weyl  expansion  of  the  Green’s  function,  and  since  the  amplitude  and  phase  of 
the  Weyl  expansion  depend  upon  the  depth  zj  —  Zj,  an  incorrect  depth1  estimate  produces 
incorrect  values  of1  the  reconstructed  optical  properties.  However,  rough  estimation  of  the 
deptli  information  can  be  tolerated  if  it  is  desirable  to  reconstruct  contrast  images  of  multiple 
objects. 

Rqi.  (2a)  reveals  h'ow  the  heterogeneity  function  and  hence  the  reconstructed  optical 
properties  of.  the  heterogeneities  vary  with  the  estimated  depth.  Choice  of  a  too  small 
depth  underestimates  the  optical  properties  and  a  too  large  depth  overestimates  the  optical 
properties.  Fig.  6(a)  shows  the  reconstructed  absorption  coefficient  of  a  spherical  object 
versus  the  estimated  deptlh  Z4  —  zaijj.  In  this  case  we  have  a  spherical  object  of  0.5  cm  radius 
2  cm  below  the  detection  plane,  i.e.,  z&  —  z0t,j=2  cm.  The  true  optical  property  variations 
of  the  spheres  with  respect,  to  the  background  are  S/ia  =  0.02  cm-1  and  8f/s  =  0.  We 
find  that  the  reconstructed  absorption  increases  as  the  estimated  object  depth1  increases. 
Tn  Fig.  6(h)  we  plot  the  ratio  of  the  reconstructed  absorption  coefficients  of  two  spherical 
absorbing  objects  versus  the  estimated  depth.  One  sphere  of  8fxa\  =  0.04  cm-1 

and  6fj,\ ,  =  0  cm"1  is  at  Q2,  1,  3)  cm  and  the  other  sphere  of  <5pa2  =  0.02  cm"1  and  8/i's2  —  0 
cm"1  is  at  3}  cm.  Two  spheres  have  the  same  size  (0.5  cm  in  radius)  and  they  are 

chosen  to  to  he  at  the  same  depth,  e.g.  2  cm  below  the  detection  plane.  Therefore  any 
depth  estimate  is  either  correct  or  incorrect  for  both  objects  at  the  same  time,  and1  we  do 
not  haw  to  take  into  account  the  additional  complexity  shown  in  Fig.  6(a).  We  find  that 
the  ratio  of  the  reconstructed  absorption  coefficients  is  not  sensitive  to  the  depth  estimation, 
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and  therefore  the  incorrect  depth:  estimate  for  contrast  image  can  he  tolerated  in  this  case. 

The  image  quality  is  also  affeoted  by  the  choice  ofi  a  pmoni,  depth.  Recall  the  heterogeneity 
function  and  the  scattered  wave  in  K-spaoe  is  coupled  to  each1  other  via  the  Weyl  expansion 
(see  Eq.  (25)).  The  factor  e~im^d~Zob^ fm  increases  exponentially  with  the  (21,/  —  21 „y).  The 
noise  (numerical  and  experimental)  can  be  amplified  at  greater  depths  (^  -  Zobj)  ■  A  series  of 
reconstructed  images  with  different  depths  are  shown  in  Fig.  7.  In  this  example  an  absorbing 
spherical  object  is  at  (2,  1,  3)  cm  and  the  scattered1  wave  is  measured  in  the  plane  at  2=3 
cm  over  a  9  x  9  cm2  square  with  steps  of  0.6  cm.  The  images  (b-f)  are  reconstructed  with 
assumptions  of  the  depth  ( zd  —  zobj)  to  be  respectively  4,  3,  2,  1,  0  cm.  We  find  that  the 
image  quality  gets  worse  (e.g.  noisier)  at  greater  depths.  The  depth-dependent  noise  and 
the  monotonia  variation  of  the  image  sharpness  make  it  difficult  to  estimate  the  true  object 
depth  from  image  sharpness.  For  a  spatially  extended  object,  however,  a  choice  of  a  shallow 
depth  is  often  sufficient  to  reconstruct  fairly  well  the  spatial  margins  of  inhomogeneities. 

Tu  order  to  obtain  better  3-D  information  with  this  diffraction  tomography  technique,  one 
cam  use  a  secondary  localization  scheme  to  deduce  the  object  depth.  An  example  would  be 
to  scam  the  phased-array  in  two  orthogonal  planes  [22,23].  Alternatively  as  shown  in  Fig.  8, 
if  we  take  two  planar  measurements  along  two  different  directions  of  the  same  sample,  the 
projection  image  1  from  the  first,  measurement  in  one  plane  (plane  1)  will  provide  the  depth 
information  for  the  projection  image  2  from  the  second  measurement  in  the  other  plane 
(plane  2). 


D.  Experimental  Results 

To  demonstrate  the  experimental  feasibility  of  this  algorithm,  we  have  performed  ampli¬ 
tude  and  phase  measurements  in  a  paralleli-plane  geometry  (Fig.  2).  We  used  a  rapid  homo- 
dyne  detection  system  based  upon  In-phase/Quadrat.ure  (IQ)  demodulation  techniques  [24]. 
The  block  diagram  of  the  experimental  setup  is  shown  in  Fig.  9.  Source  light  intensity  is 
modulated  at  100  MHz  and  the  source  power,  is  about.  ~3  mW  at  786  rim.  The  source  light 
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is  guided  into  a  large  fish  tank  of  5:0  Liters  0.75%  Intralipid  (^=0.020  era-1,  /z'=7..3  era-1) 
with  a  fiber  bundle  of  3  mm  in  diameter.  The  large  volume  of  Initralipid  enables  the  use 
of  infinite  medium  boundary,  conditions.  A  detection:  fiber  bundle  also  of  3  ram  diameter, 
couples  the  detected  diffusive  wave  to  a  fast  avalanche  photo-diode  (APD). 

The  experimental  geometry  is  shown  in  Fig.  2.  The  source  position  was  fixed  and  taken 
to  be  the  origin  of  our  coordinate  system.  As  shown  in  Fig.  2.  we  “made”  the  detection 
plane  by  scanning  a  single  detection  fiber  over  a  square  region  from  (-4.65,  -4.65,  5.0)  cm  to 
(4.65,  4.65,  5.0)  cm  in  a  plane  at  zd=5.0  cm  in  steps  of  size  Atx=Ay=0.3  cm.  The  amplitude 
and  phase  of  the  DPDW  was  recorded  at  each  position  for  a  total  of  1024  points.  Each 
data  point  takes  about  half  second.  We  directly  measured  the  amplitude  and  phase  in  the 
homoyenmus  medium  to  obtain  the  incident  wave  <I>0(!r<j). 

Tm  tliis  experiment,  a  absorbing  slice  with  dimensions  1.5  x  1.5  x  0.4  cm3  was  submerged 
in  the  turbid  medium  (0.75%  Intralipid)  at  position  (1-1.6,  -0.3,  3.0)  cm.  The  slice  was  made 
oft  resin  plus  Ti0.2  and  absorbing  dye.  Ti02  particles  (from  Sigma)  cause  the  scattering  and 
the  absorbing  dye  (900NP  from  Zeneca)  causes  the  absorption.  The  absorption  coefficient 
of.  the  slice  was  pbafibj  ~  0.20  cm-1;  its  scattering  coefficient  was  about  the  same  as  that  of 
the  background,  i.e.  ~7.3  cm-1.  The  scattered  wave  <£sc(rd)  was  obtained  by  subtracting 
the  incident  wave  $o(r.d)  from  the  measured  (total)  signal  <fi(rd). 

For  image  reconstruction,  we  first  take  the  2-D  Fourier  transform  of  the  scattered  wave 
4l^,(nci|)'  measured  at  the  detection  plane  2  =  zd.  Using  Eq.  (25)  along  with  a  priori  informa¬ 
tion  about  the  slice  depth,  we  then  obtain  the  heterogeneity  function  in  K-space  T(p,  q,  z„bj) 
in  the  plane  containing  the  slice  at  2  =  zobj .  During  this  step,  an  “m-cut”  filter  is  used  to 
neglect  high  spatial  frequency  components  with  Im{rn)  >  3.5/m(7o)  in  the  heterogeneity 
function  T.(p,  q,  zohj}.  We  then  take  2-D  inverse  Fourier  transform  of  T(p,  q,,  zobj)  with  respect 
to  spatial  frequency  (p.  q)  to  obtain  the  heterogeneity  function  T(x,  y,  zobj)  in  real  space.  Fi¬ 
nally:  we  divide  the  heterogeneity  function  T(x,y,  zobj)  by  the  background  field  $o(x,y,  zobj) 
in  the  plane  containing  the  slice  at  z =zob;i  to  obtain  a  spatial  map  of  the  reconstructed 
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absorption  variation,  e.g.,  Sfia(x,  y.  a,,y)  =  TA ^(x.  y.  ^y)/(,— $o(®>  V-.  %y))-  The  homoge¬ 
neous  background  field  $q(.t,  y,  znbj]  is  calculated,  using  the  best  estimated  optical  properties 
(/ia0=0.017  cm"1  and  /40=7.21  cm"1)  by.  fitting  the  inoidenit.  wiavie  <T>0(:rid)  measured  ini  the 
detection  plane  z =zd  to  the  exact  solution  of  DPDW’s  (Eq.  (10)).  The  reconstructed  images 
of  the  slice  are  shown  in  Pig.  10.  The  reconstructed  x-y  position  was  about  at  (-1.80.  -0.25) 
cm,  close  to  the  true  x-y  position  at  (-1.6,  -0.3)  cm.  Inaccuracies  in  the  position  measure¬ 
ments  might  account  for  the  discrepancy.  The  reconstructed  absorption  coefficient  is  well 
above  the  background  noise  level  and  close  to  the  true  value,  e.g.,  S^bj  ~  0.125  cm"1. 
Errors  in  our  estimate  of  the  slice  depth  and  background  optical  properties  estimates  are 
the  main  sources  of  inaccuracy  in  reconstructed  absorption  properties.  The  refractive  index 
mismatch  between  the  object  (~  1.46)  and  background  medium  (~  1.33)  also  contributes  to 
the  inaccuracy  in  the  reconstruction.  The  complete  reconstruction  based  upon  forward  and 
inverse  FFT.  calculations  takes  less  than  0.2  second  CPU  time  on  Sun  SparclO  workstation. 

The  feasibility  of  this  diffraction  tomography  algorithm  for  multiple  slices  and  multiple 
spatially  extended  optical  heterogeneities  has  also  been  experimentally  verified  [18].  We 
found  that:  the  optical  properties  of  thin  objects  or  the  relative  optical  properties  of  spatially 
extended  objects  (contrast  image)  with  a  priori  depth  information  can  be  reconstructed  with 
a  reasonable  accuracy. 

VI.  TWO  ISSUES  RELATED  TO  ANGULAR  SPECTRUM  ALGORITHM 

The  angular  spectrum  algorithm  provides  an  approximate  relation  between  the  hetero¬ 
geneity  function  (|or  the  heterogeneities]  and  the  scattered  wave  within  the  framework  of 
the  first,  order  Porn  approximation.  T'rt  addition  tlo  this  first  order  approximation,  it  also 
requires  the  information  ofifhe  background  optical  properties.  The  resultant  images  are  2-D 
photographic- type  images.  Tm  this  section,  we  will  consider  the  possibility  of  simultaneous 
reconstruction  of  the  absorption  and  scattering  coefficients,  and  we  will  explore  methods  to 
extract  the  background!  optical  properties  from  a  single  measurement  on  a  heterogeneous 


sample. 


A.  Absorption  and  Scattering 

So  far  we  have  assumed  that  we  have  either  purely  absorbing  inhomogeneities  or  purely 
scattering  inhomogeneities,  but  not  a  mixture.  Wie  introduce  a  dual  modulation  frequency 
approach  as  a  means  to  reconstruct  the  absorption  and-  scattering  coefficients  simultaneously. 

When  both  absorption  and  scattering  variations  are  present,  the  heterogeneity  function 
is  (see  Eiq.  (Dl)  and  Eq.  (12)) 

srr  p 

T(r)  =  <$0(r)  6pa{ r)  + - -  <Mr)  <VJr0  ■  (26) 

JJq  V 

Witih'im  a  “thin”  slice  approximation,  the  heterogeneity  function  T(r)  in  the  plane  at  z  =  zol)j 
cam  he  obtained  using  the  angular  spectrum  algorithm.  Dividing  T(j r)  by  the  incident  wave 
d>0 (]n)|  im  the  plane  at  z  =  zobj,  we  obtain  the  following  quantity,  denoted  by  F(uj),  which  is 
a  fumcfioni  of  Sjj;,a,  4 //, ,  as  well  as  the  modulation  frequency  w,  i.e. 

F&)  -  =  — 3/4  <Wr)  +  [— 3//a0  +  *^]  <Wr)  ■  (27;) 

Note  that,  the  scattering  variation  5p!s  appears  along  with  the  modulation  frequency,  while 
dfini  does  mot.  Therefore,  if  we  measure  the  scattered  wave  at  two  different  modulation 
frequemcies  a;,  amdl  cu2,  the  difference  between  the  two  the  reconstructed  F^)  and  F( u>2) 
will  only;  be  related  to 

F(ux2)  -  FM  =  i  3  —  ~  5p's  .  (28) 

Sp'g  cam  he  determined  from  Eq.  (28).  Then  by  substituting  the  resultant  Sp!s  into  Eq.  (27), 
wie  cam  them  determine  the  absorption  variation 

To  demonstrate  the  feasibility  of.  this  approach,  we  simultaneously  reconstruct  the  ab¬ 
sorption  and  scattering  coefficients  of  a  generic  slice  using  simulated  data-.  The  simulation 
geometry  is  similar  to  the  experimental  geometry  sWown  in  Fig.  2.  A  1  x  1  x  0.3  cm3  slice  of 
//„  =  0.04  cm'1  and  /i(  =  12.0  cm'1  is  placed  at  (U,  -1,3)  am.  The  source  is  at.  (0,0,0)  cm 
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and  the  homogeneous  background  has  optical  properties  ofi/i^o  =  0-02  cmr1  and  /i',0  =  8.0 
ora'1.  Note  that  the  slice  has  both  absorption  and  scattering  variations  with  respect  to  the 
homogeneous  background.  The  total  and  background  diffuse  photon  density  waves  at  zd= 
5  cm  are  calculated  for  /  =  70  MHz  and  f  =  140  MHz  using  the  finite  difference  method. 
The  reconstructed  images  are  shown  in  Fig.  11  (bj  and  (o).  The  reconstructed  absorption 
and  scattering  coefficients  are  «  0.025  cm-1  and  5p's  «  3.32  cm-1.  We  find  that  this 
approach  provides  simultaneous  estimates  of  the  absorption  and  scattering  coefficients  with 
a  reasonable  accuracy. 

B.  Extraction  of  Background  Optical  Properties 

Trniage  reconstruction  requires  knowledge  of  the  optical  properties  of  the  homogeneous 
background  medium.  For  example,  the  complex  spatial  frequency  m  =  [ko-(2n)2(p2+q2)]l/2 
in  Rq.  (25)’  depends  on  the  incident  photon  density  wave  number  k(>,  and  k0  in  turn  depends 
upon  the  absorption  and  scattering  coefficients  of  the  background  turbid  medium  (e.g., 
k0  =  |(-v/zft0  +  i  io.)/D0y/2).  Ideally  we  would  like  to  determine  the  background  optical 
properties  from  a  single  data  set  measured  on  a  heterogeneous  medium.  One  simple  way  to 
achieve  this  goal  is  to  fit  the  heterogeneous  data  set  with  a  homogeneous  model  to  obtain  bulk 
average  values  of  the  optical  properties.  We  find  that  the  results  for  this  case  are  generally 
unsatisfactory.  Fig.  12(h)  shows  the  total  photon  density  wave  T(r)  (the  amplitude,  for 
example)  from  the  absorbing  slice  experiment  where  the  detector  was  scanned  along  a  line 
symmetrically  with  respect  to  the  source.  When  fitting  all  the  data  points  with  a  simple 
homogeneous  model,  we  find  that  the  resultant  absorption  and  scattering  coefficients  are 
^*o=0.012  cmr1  and  /i'^*=C.27  cm'1,  while  the  expected  values  for  0.75%  Intralipid  are 
/irt,o=0.020  cm-1  and  /i(0=7.3f)  cm-1. 

We  can  improve  the  results  by  considering  the  symmetry  of  our  detection  scheme.  As  we 
recall  (see  Fig.  2),  our  scanning  geometry  is  mirror  symmetric  with  respect  to  the  source.  In 
Fig.  12(a),  we  project  the  3- lb  geometry  into  2-D.  to  re-emphasize  this  mirror  symmetry.  If  the 
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medium  is  homogeneous,  the  data  should  he  symmetric  with  respect  to  source;  ifithe  medium 
is  heterogeneous,  the  left-right  symmetry;  will  h'e  broken.  This  broken  symmetry  enables  us 
to  identify  the  data  points  which  are  mostly  perturbed  by  the  imbomiogeneifies.  Since  the 
phase  of!  diffuse  photon  density  waves  is  not  as  sensitive  to  the  absorption  variation  as  the 
amplitude,  it  is  safe  to  use  only  the  amplitudes  of  the  photon  density  waves  in  identifying 
the  most  perturbed  data  points  from  the  symmetry  consideration.  If  thle  left-right  difference 
in  amplitude  signals  is  greater  than  the  system  noise  level,  we  call  those  data  points  the 
most  perturbed  data  points  (See  Fig.  12(c)).  We  then  exclude  these  perturbed  data  points, 
and!  fit  the  rest  of  data  points  (both  amplitude  and  phase)  to  a  homogeneous  model  (see 
Fig.  12(d)).  We  find  that  resultant  optical  properties  are  indeed  improved,  e.g.,  =  0.015 

cm-1  and!  —  7.23  crn-1,  where  the  inaccuracy  decreases  from  ~40%  to  ~25%  in  /xa0 
and  from  118%  to  2%  in  /4o- 


VII.  SLAB  AND  SEMI-INFINITE  GEOMETRIES 


Recall  that!  the  total  photon  density  wave  $(r)  for  a  turbid  medium  with  boundaries  is 
given  by  Rq.  (14),  i.e. 


<fer'i  =  ^  /: SI(K)  G( r, if)  Sr’  +  JT(r')  G{ r,  if)  Sr! 


V' 


(29) 


On  the  surface  ofi  the  turbid  medium,  the  diffuse  photon  density  wave  satisfies  the  zero 
pflrtiai  mrrent  boundary  condition  [25]) 


4>.(n]  + 


1  +Rq£1  2 Do  r)  _  n  , 
n  bn* 


d&(r) 

dn' 


=  —  cv$(r)  , 


for  r  on  the  surface  . 


(30) 
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Here  n'  is  the  surface  normal  pointing  outward  from  the  scattering  medium,  a  —  j 
where  Refj  is  the  effective  reflection  coefficient  2.  Using  this  zero  partial  current  boundary; 
condition  in  Ecj.  (29),  we  obtain  a  general:  solution  for  the  total  diffuse  photon  density  wave 
d>(r]  in  a  finite  turbid  medium: 

<I>(r)  =  -fr-  J'  S(r')  G(r,  r')  dV  +  /T(r')  G(r.r')  dV 
0  v  v 

-/4.(r')[aG(r,rO+^]d^.  «3t) 

s 

How  the  scattered  wave  is  related  to  the  heterogeneity  function  in  this  case?  As  we 
discussed  at  the  end  of  Section  III,  the  surface  term  (the  last  term  in  the  above  equation) 
depends  on  thle  total  photon  density  wave  $(r),  and  therefore  the  surface  term  includes 
both  an  incident  wave  component  and  a  scattered  wave  component.  Analytic  separation 
of  fire  incident  wave  component  from  the  scattered  wave  component  in  the  surface  term  is 
generally;  not  feasible  though  perturbative  approaches  may  be  used  approximately. 

The  approach  wie  take  here  is  to  find  an  appropriate  Green’s  function  so  that  the  surface 
term  is  zero  by  requiring 

o:  Gfri,  r')  +  _  q  r  is  on  the  surface.  (32) 

On’ 

Note  fliat  this  boundary  condition,  as  we  discussed  in  Section  III,  is  naturally  satisfied  for 
an  infinite  turbid  medium  (no  photons  reach  the  infinity).  By  requiring  the  Green’s  function 
to  satisfy  Rq.  (32].  we  then  have  the  total!  photon  density  wave  4>(r): 

*(»]  =  ^/| A(r')  G(. r, r')  dV  +  j T(r')  G(r, r')  dV  ,  (33) 

0  v  v 

from  which  we  can  obtain  the  scattered  wave  <Rsc(|r): 


2Tli:e  exant  expression:  of  R(iff  was  derived  by  Haskell,  Tromberg  and  their  coworkers  [25] .  An 
approximate  expression  offered  by  Qnoenliuis  and  coworkers  [26],  is  in  agreement  with  the  exact 
R„ff  within  10%,  The  approximate  expression  is  R-e.ff  —  —1.440?)” 2  +  0.710??”1  +  0.668  +0.0636rt 
where  the  relative  index  of  refraction  n  =  njnituT,bid/naut,airr 
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<Mr)  =  $(rl  -  *o(»)  =  /  T('r')  GfaiO  .  (•34) 

v 

Our  task  is  to  find  the  appropriate  Green’s  function  wihli clil  satisfies  Rq.  (13)  and  the 
boundary  condition  given  by  Eq.  (32),  i.e. 

(V2  +  *(j)G(r,  r')  -  -£(r,  r')  ,  (33) 

and 

a  G( r,  r')  +  =  0  ,  for  r  on  the  boundaries.  (36) 

We  expect  the  Green’s  function  for  a  finite  medium  to  include  the  Green’s  function  in  an 
infinite  medium  G0( r,  r'),  and  an  additional  term  G/,.( r,  r')  which  results  from  backreflections 
at  the  boundaries,  i.e. 

G(r,  r')  -  G0(r,  r')  +  Gh(\ r,  r')  ,  (37) 

where  (70(jn;r'|  =  ^ .  Gfi(r,r')  is  required  to  satisfy  the  homogeneous  Helmholtz 

equation! 

(V2  +  k‘20)Gh(r,r’)=0,  (38) 


and  the  following  boundary  condition: 
dGu{ r.n') 


ot  GH( r.  n') 


d.n' 


-  -to:G0(r,r')  + 


dGo  (r,rf). 
dn'  J 


for  r  on  the  boundaries.  (39) 


A.  Slab  Geometry 

Boundaries  of  arbitrary  shapes  a, re  in  general  difficult  to  incorporate  into  the  solution  of 
the  photon  diffusion  equation  (Eq.  (14)).  Here  we  consider  a  slab  geometry  shown  in  Fig.  13. 
Within  the  slab  is  the  scattering  medium  and  outside  the  slab  is  air.  This  slab  geometry 
is  to  approximate  the  compressed  breast  configuration  which  is  suitable  for  clinical  breast 
lesion  diagnosis. 
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Suppose  the  two  surfaces  of  a  slab  turbid  medium  are  at  a  =  Zq  and  21  =  zd  as  shown  ini 
Pig.  13.  Again  we  use  the  angular  spectrum  representation  of  the  Green’s  function  Gfoflri,  ri')|. 
i.e. 

G„(r,r')  =  JfdpdqGK(p,q,z,z')  ,  (40) 

Substituting  this  equation  in  Eq.  (38),  we  find  for  any  given  spatial  frequencies  (fl,?), 
(7.wfl p,q,z,z')  satisfies  the  following  1-dimensional  homogeneous  Helmholtz  equation: 

[]kz  +  q> *0  =  0  >  (411) 


wiliere  m  =  -  (27r)2(p2  +  g2)]1/2  and  Im{m )  >  0.  The  boundary  conditions  given  by 

Fkj.  (39)  for  a  slab  geometry  shown  in  Fig.  13  can  be  rewritten  for  the  angular  spectrum 
fifcCPs  ?:asz*)  as 


aCH{p:  4,  z,  z'  =  z0)  -  - £°1  =  —(a  +  im)G0{p,  q,  z,  z'  =  z0)  , 

a(7/i(p.  q.  z.  zf  =  2(f)  -H  — h^’ ^LJ. - —  =  —(a  -H  irn)Go(p,q,  z,  z'  =  2,/)  , 

where  G^q.z.z')  =  gjven  by  the  Weyl  expansion  (see  Eq.  (19)). 

The  general  solution  of  Gh{p,  q,  z ,  z')  has  the  form  of 


(42) 

(43) 


G,((Mjz)z')=.4e-'  +  Br 


(44) 


The  first  term  represents  the  wave  which  is  reflected  by  the  lower  surface  at  2  =  z0  and  then 


propagates  forward,  along  +z  direction,  i.e.,  the  “transmission”  component;  the  second  term 


represents  the  wave  which:  is  reflected  by,  the  upper  surface  at  z  —  zd  and  then  propagates 
backward  along  -21  direction,  i.e..  the  “reflection”  component.  Coefficients  A  and  B  can  then 
be  solved  using  the  boundary  conditions  given  by  Eqs.  (42)  and  (43).  After  some  algebra, 
we  find  that. 


A\  =  haimx  +  he' 


B  =  heimz  +  he 


(45) 


where  f\,h,h  and  /*  are  given  by, 
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wiitlhl 


}\  =  4(«2  +  ™2)  e-im(Zrf+2o);  /2  =  -4(«  +  imf  dm[zd~z^  ,  (4.6) 

u  P. 

f9  = -k  (a  +  imf  eiM*d-*> ) ,  f«  =  j  ('a2  +  mi2)  eim(2d+2o)  ,  (47) 

/■„  =  —  ,  P  =  (a  +  imfeim(zd~za)  -  (a  -  xmJze"im(^"aDj  .  (48) 

2m 


Pm  ally  for  a  slab  geometry,  the  Fourier  component  of  total  Green’s  function  Gl(p,  <?,  a,  a*) 


im  K -space  is 


G(p,  9,  z,  a/)  =  G0(p, z')  +  Gh(p„  q,  z ,  z')  = 

yAm\z-z’\  +  jieim{z+z')  +  fa-imt.z-z')  +  f^m{z-z')  +  } 


(49) 


Using  G{p,q,z,z')„  the  relation  between  the  scattered  wave  <3>sc(p,  q,  zd)  and  the  hetero¬ 
geneity;  function  T(p,  q,  zobj)  for  a  slab  geometry  within  “thin”  slice  approximation  is 

rfu  S  ®3cfaq,Zd)  ^ 

Tip,  q,  zobj)  «  — -s- - r  ■  (&0J 

AzG{p,q,  zob,j) 

The  2-D  inverse  Fourier  transform  of  T(p,  q,  zobj)  gives  the  heterogeneity  function  T(x,  y,  zobj) 
in  real  x-yi  space.  The  optical  properties  of  the  inhomogeneities  can  then  be  obtained,  e.g., 
for  absorbing  objects,  we  have  Spa^x., yy  z^j)  =  Tabs{x,y,zobj)/i—^  ^ab{:x,  y,  zobj)):  and  for 
scattering  objects,  wie  have  dlp/Jr,?/,^)  =  Tsc{x,  y,  zabj)/{ 3D° ^  %lab(x,y,  zobj)). 

Using  the  appropriate  Green’s  function  (Eq.  (49))  for  a  slab  geometry,  we  reconstructed 
a  2-D  optical  image  of  a  slice  embedded  in  a  slab  turbid  medium.  The  slab  geometry  is 
shown  in  Fig.  13  where  the  two  surfaces  are  at  planes  z  =  0  and  z.  =  5  cm,  respectively. 
The  source  is  at  the  origin.  A  1.0  x  1.0  cm2  slice  of  0.3  cm  thick  is  at  (1,-1, 3)  cm  (see 
Fig.  14  (a))  and  the  detection  plane  is  at  the  top  surface  of  the  slab(  z  =  5  cm).  The 
slice  has  a  higher  absorption  coefficient,  than  the  background  medium  but  shares  the  same 
scattering  coefficient:  with  the  background,  e.g.,  pfaobj  =  0.04  cm-1  and  p!sobj  =  8.0  cm  1 
for  the  slice  and  n„ 0  =  0.02  cm-1  and  p\'sCl  =  8.0  cm-1  for  the  background.  The  total  and 
background  diffuse  photon  density  waves  at,  the  top  surface  z=5  cm  are  calculated  using  the 
finite  difference  method. 
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The  reconstructed!  absorption:  image  using  th'e  appropriate  Green’s  function  (jEq.  (]49]|]| 
for  the  slab  geometry  is  shown  in  Fig.  14'  (b)  .  and  the  reconstructed  absorption  is  « 

0.0240  cm-1,  which  is  close  to  the  expected  value  Sgf  =  b.0200  cm*1.  For  comparison, 
we  also  reconstructed  the  image  of  the  slice  using  thle  wrong  Green’s  function,  i.e.  the 
Green’s  G'0(r, r')  which  is  only  correct  for  the  infinite  medium  (Eq;.  (|19j).  The  resultant 
absorption  image  is  shown  in  Fig.  14  (c).  We  found  that  the  position  of  the  slice  can  be  well 
reconstructed  by  using  different  Green’s  functions.  However  the  image  shown  in  Fig.  14  (c), 
which  is  reconstructed  by  using  the  wrong,  infinite  Green’s  function,  has  more  artifacts  than 
that,  shown  in  Fig.  14  (b),  which  is  reconstructed  by  using  the  right  slab  Green’s  function. 
Notice  that  the  Green’s  function  for  a  slab  geometry  \G(p,  q,  z,  z')  \  is  smaller  than  the  Green’s 
function  for  an  infinite  geometry  \G0{p,q,  z,  z')\  simply  because  we  lose  photons  through  the 
finite  boundaries.  When  we  use  the  Green’s  function  of  an  infinite  medium  to  reconstruct 
the  image  for  a  slab  geometry,  the  overestimate  of  the  Green’s  function  is  responsible  for  the 
nosier  image  structures  (artifacts)  in  Fig.  14(c).  The  overestimated  infinite  Green’s  function 
also  results  in  smaller  reconstructed  optical  properties,  e.g.,  the  reconstructed  value  by  using 
the  wrong  infinite  Green’s  function,  8pTa“rong  =  0.0056  cm-1,  is  about  4  times  as  small  as  the 
value  6ff™ioh  =  0.0240  cm*1  reconstructed  by  using  the  appropriate  slab  Green’s  function. 
We  see  that  the  appropriate  Green’s  function  for  a  slab  geometry  (Eq.  (49))  produces  cleaner 
images  and  more  accurate  optical  properties  than  the  Green’s  function  which  is  only  suitable 
for  an  infinite  medium  (jEcj.  (119)). 

B.  Semi-infinite  Geometry 


As  an  extension  of  the  above  derivation,  we  can  easily  obtain  the  Fourier  component  of 
the  total  Green’s  function  0™mi(p,  q,  zi,  z1)  for  a  semi-infinite  turbid  medium.  Alternatively 
we  can  start:  with  the  Green’s  function  for  a  slab  geometry  (  Eq.  (49)),  then  move  the 
lower  boundary  of  the  slab  in  Fig.  13  to  the  negative  infinity,  i.e.  zq  -4  — co.  Note  that 
Intern)  >  0  and  therefore  all  terms  in:  Eq,  (49)  with  rr’"1*0  vanish  when  z0  -4  -oo.  The 


Fourier  component  of  the  total  Green’s  function  d?MM*(|p.  q,  a,  21']  for  a  semi-infinite  medium 
at  any  spatial  frequencies  (p,  q'j  in  K-space  is  thus 

Gsemi(p,  q,  2,  z1)  =  f0  eimlz~^  +  e"im(x+2,)  .  (:31) 


Here  the  coefficients  /o  and  /|erm  are  given  by 

fo  =  fi,  tr“  =  -/»  e*™"  .  (52) 

2  m  a  —  vm 

The  first  term  on  the  right  hand!  side  of  Eq.  (51)  represents  the  Green’s  function  in  an  infinite 
medium.,  and  the  second  term  represents  the  wave  which  is  reflected  by  the  boundary  at 
zi  =  a,/  and  propagates  backward  along  the  negative  2-direction.  For  arbitrary  boundaries 
the  solution  of  the  Green’s  function  are  in  general  difficult  to  obtain. 


C.  Re-emission  Geometry 

Tn  the  previous  discussions  the  source  and  the  detector  were  assumed  to  be  on  the 
opposite  sides  of  the  inhomogeneity.  This  configuration  is  called  transmission  (see  Fig.  15 
(]a||.  Tt.  is  suitable  far  two-plate  soft  compression  geometry  in  breast  cancer  studies,  with  the 
source  placed  on  one  plate  and  the  detector  scanned  over  the  other  plate.  Interestingly  the 
derivation  is  not  limited  to  this  transmission  configuration.  Recall  that  dependence  of  the 
angular  spectrum  algorithm  on  the  source  position  is  implied  in  the  heterogeneity  function 
(see  Eqs.  (11  )|  and  (12)|j|.  The  relation  between  the  heterogeneity  function  and  the  scattered 
wave  measured  at  the  detection  plane  (see  Eq.  (25))  does  not  explicitly  depend  on  the  source 
position.  The  light  source  and  the  detector  can  be  placed  on  the  same  side  or  on  the  opposite 
sides  of  the  object  without  affecting  the  conclusion  of  the  above  derivation.  Hence,  we  can 
apply  the  algorithm  equally  well  to  another  geometry  -  th'e  re-emission  geometry  (see  Fig.  15 
(hJJ  [|27!.28].  Tn  the  re-emission  configuration,  the  detector  scans  in  the  plane  which  contains 
the  source.  This  could  he  necessary,  for  exam  pile,  in  brain  function  studies.  The  re-emission 
geometry  could  also  he  useful  for  studies  of  large  dense  breast  tissues  in  which  fewer  photons 
pass  through  the  tissue. 


In  the  transmission  geometry,  we  measure  the  scattered'  wavie  propagating  forward,  away 
from  the  source;  in  the  re-emission  geometry,  we  measure  the  scattered  wavie  propagating 
backward  towards  the  source.  For  a  re-emission  geometry  and  within  a  thin  slice  approxi¬ 
mation,  the  relation  of  the  heterogeneity  function  f(p,  qn  zoh:j)  in  K-space  with’  the  measured 
scattered  wave  in  the  plane  z  —  zd  is  given  by  the  same  equation  as  for  the  transmission 
geometry  (Eq.  (25)).  Here  we  rewrite  the  relation  for  the  re-emission  geometry: 


T(p,  q,  z0bj)  = 


$»c(p,q,2d) 

A z  G(p,  q ,  zd,  zobj) 


where  an  appropriate  Green’s  function  for  an  infinite  medium  (Eq.  (19))  or  a  slab  medium 
(Rq.  (4.9))  has  been  assumed. 

Simulations  hla/ve  shown  the  applicability  of  the  algorithm  to  the  re-emission  geometry. 
Consider  an  absorbing  spherical  inhomogeneity  of  0.5  cm  radius  at  (2,  1,  2)  cm  embedded 
in  an  otherwise  homogeneous  slab  turbid  medium  (see  Fig.  15  (c)).  The  two  surfaces  of 
the  slab  are  at  za  =  0  cm  and  zd  =  4.0  cm,  respectively.  The  absorption  and  scattering 
coefficients  of)  the  sphere  are  p,,— 0.04  c.m_1  and  /i'.=8.0  cm-1  while  the  background  optical 
properties  are  /-Wi^O-02  cm"1  and  p'o=8.0  cm-1.  For  the  re-emission  configuration,  both 
the  source  and  detector,  are  placed  on  the  top  surface  of  the  slab,  i.e.  in  a  plane  at  zd  =  4 
cm.  The  scattered  wave  is  calculated  over  a  9  x  9  cm2  region  with  x-y  steps  of  0.6  cm.  The 
source  is  placed  at.  the  center  of  the  square  scanning  region  in  the  detection  plane,  i.e.  at  (0, 
0.  4')  cm.  The  reconstructed  image  for  the  re-emission  geometry  is  shown  in  Fig.  15  (d).  For 
comparison  we  also  reconstruct  the  image  of  the  same  object  for  the  transmission  geometry. 
Tn  this  case  the  source  is  at  the  origin  (0.  0,  0)  c.m,  on  the  lower  surface  of  the  slab,  with 
all  other  configurations  kept  the  same  as  in  the  re-emission  geometry.  The  image  is  shown 
in  Fig.  15  (]e|.  We  see  that,  the  image  quality  in  these  two  configurations  is  about  the  same. 
The  ratio  ofitlie  reconstructed  absorption  coefficient,  for  the  re-emission  geometry  to  that  for 
t.be  transmission  geometry  is  Spr^_em/6f/a^ans  ~  IU .  The  finite  object  size  (as  opposed  to 
a  “tliin”  slice)  might  contribute  to  the  small  difference  in  the  reconstructed  absorption. 
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VIII.  SUMMARY 


We  have  presented  a  full  exposition  of  our  recent  work  that,  employs  the  angular  spectrum 
algorithm  for  optical  diffraction  tomography  with  diffuse  photon  density  waves.  The  image 
reconstruction  becomes  practically  easy  for  thin  heterogeneities  wherein  the  heterogeneity 
function  of!  interest  is  proportional  to  the  scattered  wave  measured  at  tire  detection  plane, 
i.e.  T(p,  q,  z0bj)  oc  <hsc(p,  q,zd).  We  have  shown  that'  although:  this  relation  is  accurate  only 
for  thin  in  homogeneities,  it  provides  an  approximate  short  cut  for  fast,  2-D  projection  imag¬ 
ing  of  spatially,  extended  objects.  The  reconstruction  is  very  rapid,  requiring  only  a  forward 
and  inverse  Fourier  transform,  e.g.,  it  takes  less  than  0.2  second  on  a  SparclO  workstation  to 
reconstruct  an  image  of  ~1000  pixels.  For  spatially  extended  objects,  although  the  recon¬ 
structed  optical  properties  are  not  accurate,  the  ratio  of  the  reconstructed  optical  properties 
of,  multiple  objects  are  close  to  the  true  ratio.  In  this  sense  we  say  that  contrast  image 
can  still  be  obtained  by  using  this  algorithm.  The  feasibility  for  using  this  algorithm  for 
imiage  reconstruction  of  absorbing-  and  scattering  iruhomogeneities  has  been  experimentally 
demonstrated  {18]. 

We  have  extended  the  theory  to  other  geometries  including  the  slab  and  the  semi-infinite 
geometry  for  both  transmission  and  re-emission  configuration.  The  theory  was  confirmed  in 
si  miul  at  ion  experi  miemts . 
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APPENDIX  A:  CONVENTIONS  USED  IN  THIS  PAPER  REGARDING 


FOURIER  TRANSFORM 

The  conventions  regarding  the  forward  and  inverse  Fourier  transforms  are  as  follow. 
Consider  a  function  f(x)  in  1-dimension: 

Forward  Fourier  transform: 

F(p)  =  J  f(x)  ea"”  dx  ;  (Al) 

Inverse  Fourier  transform: 

f{x)  =  I  F(q)  eriH  xp  dp  ;  (A2) 

and  the  ^-function  is  therefore  given  by 

&{p)  =  J  ei27!xp  dx  ,  h(x)  =  f  e~iHvx  dp  .  (A3) 

Using  these  conventions,  we  eliminate  the  27T  factor  outside  the  integral  of  forward  and 
inverse  Fourier  transforms. 

APPENDIX  B:  THE  WEYL  EXPANSION  OF  GREEN’S  FUNCTION 

Consider  an  infinite  turbid  medium.  paQ  and  /i(0  are  respectively  the  absorption  and 
scattering  coefficients.  The  Green’s  function  £?o(r.,r')  satisfies  the  following  equation: 

(V2  +  Ml)  G0(r,  r')  =  -<?( r,  r')  ,  (Bl) 

where  /,-.0  =  [(—'*74.0  +  i  a>)/D0]1/l2  with  Tm(kn)  >  0;  D0  =  v/ (3  is  the  photon  diffusion 
coefficient..  The  solution  of  the  Green’s  function  is  [19] 

r>fc0|r— tV| 

no(r.r')  =  — - (B2) 

4-7r|r  —  r  | 

The  Green’s  fumctioni  is  related  to  its  Fourier  transform  by 


Ot>  W  =  ///  &(p.  1:  n)  dp.  <k,  dn  , 


(B3) 


where  we  assume  r*  =  0  without'  losing  generality  and  (jp,  q,  n]  are  the  spatial  frequen¬ 
cies.  Plugging  Eq.  (B3)  into  Eq.  (1R2)  and  using  the  integral  expression  for  the  ^-function 
(Eq.  (A3)),  we  have 

JJf  G0(p,  q,  n)  [kl  -  (2tt  f(p2  +  q2  +  n?)]e~i2*  ^+®+n2)  dp  dq  dn 

—  —  Jff  e~l2w  (px+w+uz)  dp  dq  dn  .  (B4) 


Without!  a:  rigorous  proof,  we  can  obtain  the  Fourier  transform  of  the  Green’s  function  just 
by  looking  at  both  sides  of  the  above  equation,  i.e. 

1  li 


Go(p,  <?>  n)  ^27r)2(p2  +  q1  -fi  n2)  _  kfi  (2n)2n2  -  m2 
where  m  ~  [kl  —  (2ir)'2(p2  +  q2)]'12  and  Im(m)  >  0.  Eq.  (B3)  can  then  written  as 


flB5j 


i2n  nz 


(B6) 


2n)2n2  —  m2 

The  integral  over  spatial  frequency  n  can  first  be  done  by  “pole”  structure  analysis.  There 
are  two  poles  in  the  integral  over  n  as  shown  in  Fig.  16.  For  z  >  0,  we  require  Im(n)  <  0 
to  ensure  the  convergence  of  the  integral  over  n.  Therefore  we  choose  the  pole  in  the  lower 
lialfi  space  (]Fig.  16(a)),  e.g.,  n  -  -m/2 7r  (recall  Im(m )  >  0  which  gives  Im(n )  <  0).  Note 
that  the  integral  is  along  the  clockwise  direction  which  gives  us  an  extra  minus  sign.  The 
resultant  integral  is 

t.2t  nz 


f  (|2tc 


sj  Z  2,TT\T\jZ\  »(j  n 

dn  ~  —2w  i _ - _  n=-m/2^  J_  imz 

an—  z/t  l  /n  %n/  9m' 


(B7) 


W>+g:)Cn-«) . . 

Similarly,  for  a  <  0  we  require  Tm[v)  >  0.  Therefore  we  choose  the  pole  is  in  the  upper  half 
space  (] F i g .  16(]h|)|,  e.  g..  n  =  m/2?r.  The  integral  is  along  the  counter-clockwise  direction  so 
there  is  no  extra  minus  sign  in  this  case.  The  resultant  integral  is  thus 

f  g— nz  e-i2nnz 

J 


_  n=mj 2ir  _imz 

(2y]2(n+§)  2m 


(B8) 


2x>  V"  1  2t r 

Gomihining  Rqs.  (B7:|  and  (]BS]|,  we  have  the  general  expression  for  the  integral  over  n: 

y  —  i^TT  Tl 


I  tl-nV-n2  -  m2  2 rn 


(B9) 
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Substituting  this  equation  into  Rq.  (B6).  we  them  emd  up  with1  the  Weyl  expansion  of  the 
Green’s  funotion: 

Qo(r)  =  [f  d®  dq  er^ >  (BIO) 

where  rn  -  [fcjj  -  (2n )2(p‘2  +  </2)]1/2  and  Im(m )  >  0. 

The  Weyl  expansion  represents  the  superposition  of  elementary,  harmonic  waves  in  the 
x-  and  ^-directions  (e-*Mpz+w));  the  harmonic  waves  exponentially  attenuate  in  the  2- 
dlirection  away  from  the  plane  2  =  0  which  contains  the  source.  The  harmonic  waves  and 
the  attenuation  factor  hgll  are  so  combined  that  the  double  integral  in  Eq.  (BIO)  over  all 
the  spatial  frequencies  (p,  q)  yields  the  elementary  damped  spherical  wave  on  the  left-hand 
side  of  Rq.  (B10).  i.e.,  G0(r)  =  eifc°r/;(47r  r). 
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FIGURES 


Incident 


FIG.  1.  Im  the  presence  of  optical  inhomogeneities.  the  spherical  wave  front  of  the  incident 
wave  is  distorted.  The  total  photon  density  wave  is  the  sum  of  the  incident  wave  and  the  scattered 


wave. 


FIG.  2.  Illustration  of  2-D  geometry  which  we  consider  for  the  image  reconstruction  algorithm 
based  upon  K-space  spectrum  analysis.  The  scattered  wave  $sc{x>V>z<i)  (lor  Its  spatial  Fourier 
component  <i>sr(p.,  q,  z<j.)  Is  determined  at  the  detection  plane  z  —  by  scanning  the  detector 
over  a  square  region.  Without  losing  generality  we  assume  the  optical  heterogeneities  are  located 
bielowi  tlile  detection  plane  at  2  =  2 d.  A  point  source  can  be  placed  anywhere  in  the  turbid 
medium.  In  practice  the  point  source  and  the  detection  plane  are  either  on  the  opposite  side  of 
tlile  heterogeneities  (transmission)!  or  both  on  the  same  side  of  the  heterogeneities  (re-emission). 
In  this  figure  the  point  source  happens  to  be  placed  at  the  origin  of  our  coordinate  system  for 
demonstration,  of  a  transmission  measurement  geometry, 
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PIG.  3.  TUiie  heterogeneity  functions  Ti  and  T2  with  spatial  frequency  p  from  two  slices  propa¬ 
gate  to  tlile  detection  plane  at.  z~zd  where  they  add  up  to  make  the  scattered  wave  4>sc  in  K-space 
at  the  sarnie  spatial  frequency,  p. 


PIG.  4.  (a)  and  (bt)  respectively  show  the  amplitude  attenuation  and  phase  shift  associated 


Wiithi  the  Weyl  expansion  in  Kl-space  versus  spatial  frequencies  (p,q).  Note  in  (a)  the  z-axis  is  the 
log  of  tl ile  amplitude  of  a*H zd~zj)/m[  hv,  (]b)  the  z-axis  is  the  phase  of  eim^Zd~z^  f\ m  in  degrees,  (c) 
and  (d)  alihwi  tlile  amplitude  attenuation  and  phase  shift  versus  the  depth  zd  -  Zj  for  given  spatial 
frequencies  (0.1,  0.1)  cm-1  (solid  lines)  andl  (0.5,  0.5)  cm  *  (dashed  lines). 
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T(p,q,zobj) 


-  Usc(p,q,zd) 

Detection 

—  Plane 


Z=0  z  obj  Z=zd 

FIG.  5.  The  liei. exogeneities  are  considered  to  be  thin  which  locate  within  a  thin  slice  at  2  =  zol)J 
in:  parallel  to  the  detection  plane.  The  heterogeneity  function  within  this  th'in  slice  is  approximately 
uniform  and  the  heterogeneity  function  is  zero  elsewhere. 
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BIG.  6.  (a)  shows  IJlue  reconstructed  absorption  coefficient  versus  the  depth  estimation.  The 
data  points  in  (a)  are  normalized  by  the  absorption  reconstructed  at  the  depth  where  the  object 
is.  e.g.  z,i  —  znbj  ^  2  ora.  (|b!j  shows  the  ratio  of  reconstructed  absorption  of  two  spherical  objects 
versus  the  depth  estimation.  Although  the  ratio  is  only  approximately  reconstructed  (e.g.:  the  true 
ratio  is  2).  the  ratio  is  relatively  insensitive  to  the  depth  estimate. 
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FIG.  7.  IMlustratliom  of;  the  dependence  of  reconstructed  images  on  the  estimated  depth.  The 
detection1  plane  is  at  z—5  cm  and  an  absorbing  object  shown  in  (a)  is  at  (2.  1.  3)  cm  which  is  2  cm 
below  the  detection  plane,  (b-f)  are  the  images  reconstructed  with  an  estimated  depth  respectively 
at  4  omi.  3  cm.  2  cm,  1  crn  and  0  cm. 
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EIG.  8.  Illustration  of  how  Bo  obtain  a  3-D  image  from  two  projection  images  reconstructed 


from  two  measurements  along  two  orthogonal  directions.  Image  I  from  the  measurement  in  plane 
I  provides  the  depth:  inform  at  Jon  for  image  2  from  the  measurement  in  plane  2. 
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Detection  plane  1 


PIG.  9.  Block  diagram  of  the  homodyne  setup.  APD-Avalanche  Photo  Diode;  A- Amplifier; 
SSB  I/, Q- Single  Side  Band  In-phase/Quadrature-phase  Demodulator;  LO-Local  Port;  KF-Ractio 
Frequency  Port;  A  DC- Analog  to  Digital  Converter;  GPIB-General  Purpose  Interface  Bus. 
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9.30  cm 


(a)  (b)  (c) 


FIG.  10.  (a)  shows  the  exact  x-y  position  of  a  thin  absorbing  slice.  (,b)  shows  the  surface  plot  of 
the  reconstructed  absorption  variation  (S/j,Taec)  using  the  angular  spectrum  algorithm,  (e)  illustrates 
the  reconstructed  2-D  phbtograpMc  image  of  this  slice.  Agreement  between  the  reconstructed 
position  amid  the  exact  position  as  shown  in  (a)  can  be  readily  found. 
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Thin  Slice 


Absorption  Image 


Scattering  Image 


(a)  (b)  (c) 

FIG.  ID.  (a)  shows  a  thin  slice  object  at  (1,  -1,  3}  cm.  The  slice  is  0.3  cm  thick  with  its  1 
x  1  cmi2  surface  in  parallel  to  the  detection  plane  at  =  5  cm.  The  scattered  waves  at  two 
modulation  frequencies  (70  MHz  and  140  MHz)  in  the  detection  plane  at  zj  =  5  cm  are  calculated 
using  finite,  difference  method  over  a  9.3  x  9.3  cm2  region  with  x-y  steps  of  0.3  cm.  (b)  and  (c) 
show  the  absorption1  and  scattering  images  reconstructed  simultaneously  using  the  dual  modulation 
frequency  approach.  The  reconstructed  position  of  the  slice  is  close  to  its  true  position  and  the 
reconstructed  absorption'  and  scattering  properties  are  close  to  their  true  values.  See  Section  VI(A) 
for  details. 
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PIG.  12.  (a)  si ilows  a  2-D  version  of.  the  experimental  geometry  in  Pig.  2.  The  detector  scans 
along  a  line  from  left  to  right  symmetrically  with  respect  to  the  source,  (ib)  shows  the  raw  data 
measured  on  a  heterogeneous  medium  by  scanning  the  detector  along  a  line  from  left  to  right,  (c) 
shows  the  most  perturbed  data  points  far  which  the  left,  and  right  differences  are  greater  than  the 
noise  level  of  our  detection  system  (je.g.,  2.5  mV  in  this  case),  (d)  show  the  rest  data  points  after 
the  most  perturbed  data  points  are  filtered  out.  Thte  background  optical  properties  can  then  be 
obtained  by  fitting  the  data  points  shown  in  (d)  to  a  homogeneous  model. 
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FIG.  13.  A  slab  geometry  is  considered  for  the  boundary  problem.  The  slab  is  infinite  long 
hut  has  a  finite  thickness,  e.g.  —  z<j.  One  surface  of  the  slab  is  at  plane  z  —  zq  and  another 
surface  is  at  plane  z  —  zci .  The  turbid  medium  is  between  these  two  planes  and  outside  the  slab 
is  Tiou-scatteriimg  media  such  as  air.  This  slab  geometry  is  quite  suitable  for  a  compressed  breast 
configuration  in  clinical  studies. 
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FTG.  14.  (a)  shows  the  position  of  a  0.3  orn  thick  1.0  x  1.0  cm3  absorbing  slice  at  (1,  -1,  3) 
omi  in  a  slab  turbid  medium.  The  two  surfaces  of  the  slab  are  respectively  at  planes  z=>0  cm  and 
73=5  omi.  Thte  source  is  at  origin  at  one  of  the  slab  surface  (z=0  cm)  and  the  detector  scans  at  thle 
other  surface  (a  =  5  cm).  The  reconstructed  absorption  image  using  the  “slab”  Green’s  function 
(Eq.  (40))  is  shown  in  (b).  The  reconstructed  absorption  image  using  the  wrong  “infinite”  Green’s 
function  (Eq.  (19))  is  shown  in  (c). 
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Transmission 


Re-emission 


Absorbing  Sphere  Re-emission  Image  Transmission  Image 

(C)  (d)  (e) 

PIG.  15.  (la)  shows  the  transmission  geometry.  The  source  is  at  the  origin  and  the  detector 
soan«  ini  a  plane  at  z  —  z,j.  (|b)  shows  the  re-eirhssion  geometry  where  the  source  is  at  the  center  of! 
the  detection  plane  at  z  =  (]c)  illustrates  a  spherical  absorbing  object  at  (2,  1,  2)  cm  embedded 

within  a  slab  turbid  medium.  The  two  surfaces  of  the  slab  are  at  z0  —  0  cm  and  zd  —  4.0  cm, 
respectively.  For  both  transmission  and  re-emission  geometries,  the  scattered  waves  in  the  detection 
plane  at  z,i  —  4  cm  are  calculated  using  the  exact  DPDW  solution  for  a  slab  geometry  over  a  9  x  9 
cm;2  region  with,  x-y  steps  ofO.fi  cmi.  The  reconstructed  images  for  the  re-emission  and  transmission 
geometries  are  shown  in  (dj  and  (e),  respectively.  The  two  images  look  similar.  We  also  found  that 
tine  reconstructed  absorption  coefficients  are  also  about  the  same  under  both  geometries. 
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(a)  (b) 


FIG.  16.  There  are  two  poles  for  the  integral  over  n  in  E'q.  (B6).  (ja):  for  z  >0  the  singularity! 
is  at  n  —  - mj 2k  and  the  integral  is  done  along  the  lower  close  curve;  (jb):  for  z  <0  the  singularity 
is  at  n  =>  m/2n  and  the  integral  is  done  along  the  upper  close  curve. 
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Abstract: 

A  Clinical,  compression-plate  diffuse  optical  imager  is  used  to  determine  the  bulk  optical 

properties  and  blood  saturation  in  breast  tissue  of  volunteers  and  in  tissue  phantoms. 

©  1999  Optical  Society  of  America 

OCIS  codes:  (170.3830)  Mammography;  (170.5270)  Diffuse  Photon  Density  Waves;  (170.3660)  Light  Propaga¬ 
tion  in  tissues;  (170.5280)  Photon  Migration 

Diffuse  optical  tomography  (DOT)  in  the  near-infrared  (NIR)  is  currently  emerging  as  a  viable  means  for  breast 
tumor  imaging  and  specification [1].  This  method  relies  on  optical  properties  and  biological  factors  such  as  the  blood 
oxygen  saturation  to  enhance  tumor  specificity  and  sensitivity.  It  is  important  to  study  the  dispersion  of  normal 
breast  bulk  optical  properties  in  order  to  asses  the  contrast  expected  from  tumor  tissue.  For  example,  various  factors 
such  as  age  and  the  menstural  cycle  are  expected  to  influence  the  optical  properties [2]  and  should  be  characterized. 

In  this  study  we  employ  a  clinical,  compression-plate  optical  breast  imager  for  measurements.  The  instrument  uses 
three  wavelengths  -  830nm,  786nm  and  750nm,  and  employs  a  scanning,  fiber-coupled  PMT  detector  for  detection. 
The  lasers  are  coupled  to  sixteen  different  source  positions  and  a  dicon  multiplexer  is  used  to  select  the  source  position. 
In  our  clinical  measurements  so  far  we  have  been  using  single  source  position  due  to  time  constraints.  The  lasers 
are  modulated  at  140MHz  to  produce  a  diffuse  photon  density  wave  (DPDW)  in  the  medium.  The  patient  lies  in 
the  prone  position  and  her  breast  is  inserted  into  a  small  tank  filled  with  a  matching  solution  of  Intralipid  through 
an  openning  on  the  bed.  The  detector  scans  along  the  output  plate  glass  surface,  and  the  source  is  attached  to  a 
compression  plate  which  applies  a  gentle  compression  to  the  breast.  Usually  the  range  of  compression  is  4.5cm  to 
7cm.  It  takes  «  15  minutes  to  acquire  data  from  a  9.6cm  (x)  by  4.8cm  (y)  scan  region  with  153  (17x9)  points.  We 
take  two  sets  of  data  for  each  patient;  1)  tank  filled  with  Intralipid  only  ,  i.e  without  the  breast,  which  allows  us 
normalize  for  instrument  response,  2)  tank  filled  with  Intralipid  and  the  patient  breast.  For  the  last  six  months  we 
have  been  acquiring  data  at  the  Hospital  of  University  of  Pennsylvania  from  volunteers  with  normal  breasts. 

In  order  to  test  the  feasibility  of  our  instrument  and  approach  we  also  performed  phantom  measurements  using 
balloons  filled  with  different  concentrations  of  ink  and  Intralipid.  The  balloons  were  inserted  into  the  Intralipid  tank 
to  fill  a  volume  which  is  similar  to  the  average  breast  volume  obtained  from  the  clinical  trials.  A  large  fraction  of  the 
balloon  was  kept  above  the  Intralipid  level  to  simulate  the  chest  wall.  For  analysis  we  follow  the  same  procedure  as 
done  with  the  healthy  patients.  We  have  tested  nine  phantoms  with  different  optical  properties,  different  background 
optical  properties  and  different  volumes.  Here  we  show  the  results  from  786nm.  The  results  from  the  other  wavelengths 
are  comparable. 

To  estimate  bulk  optical  properties  we  solve  the  forward  diffusion  problem  numerically  using  a  finite  difference  solver 
for  arbitrary  heterogeneities  in  the  domain  [3,  4].  To  obtain  the  bulk  properties  we  specialize  this  approach  ,  segmenting 
the  image  volume  into  Intralipid  background  and  breast  volumes  (see  fig.(l)).  During  the  measurement  we  obtain  an 
outline  of  the  breast  in  the  x-y  plane,  and  the  y-z  plane  in  order  to  define  the  segmentation  geometry.  The  top  x-z 
segmentation  plane  wall,  however,  is  assumed  to  extend  far  above  the  Intralipid  level  in  the  tank  and  thus  models 
the  rest  of  the  chest.  Therefore,  we  end  up  with  a  three  dimensional  T-shaped  heterogeneity  that  approximates  the 
breast  volume  as  illustrated  in  fig.  (1)  .  The  relevant  boundaries  are  also  shown. 

An  inverse  solver  based  on  this  solution  which  in  its  simplified  version  can  be  used  to  obtain  the  bulk  properties 
of  the  heterogeneity  volume  was  developed.  First  we  divide  the  heterogeneous  data  by  the  background  data  giving 
a  “normalised”  DPDW.  The  solver  is  then  initialized  with  the  background  properties,  convergence  is  tested  as  the 
squared  difference  between  the  measured  and  calculated  normalized  DPDWs  .  If  it  did  not  converge  ,  then  the  optical 
properties  are  updated  according  to  a  rule  based  on  distorted  born  iterative  method  and  the  forward  solution  is 
tested  for  convergence  once  more.  The  process  is  repeated  until  either  a  particular  number  of  iterations  is  reached  or 
a  certain  criterion  in  the  error  is  achieved.  Here  we  simplified  this  process  for  the  bulk  properties  by  summing  up  the 
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weight  matrices  (essentially  averaging)  in  the  solver  for  the  heterogeneity  region  and  by  not  allowing  the  properties 
of  the  homogenenous  background  region  to  change  (as  opposed  to  calculating  every  single  voxel  in  the  domain).  We 
allowed  for  a  maximum  of  ten  iterations  and  stopped  earlier  if  the  error  falls  below  10  .  We  find  that  this  criterion 

is  generally  adequate  for  our  numerical  simulations  and  phantom  tests. 

To  illustrate  our  results,  we  show  four  sets  of  data  from  balloon  phantoms  as  follows  where  the  balloon  volume  and 
other  parameters  were  kept  constant  but  the  optical  properties  were  changed.  The  background  properties  were  0.05 
and  8  cm-1  for  the  absorption  and  scattering  (^coefficients  respectively. 


Expected  (cm-1) 
Absorption 

Scattering 

Calculated 

Absorption 

(cm-1) 

Scattering 

0.05 

12 

0.0587 

13.54 

0.05 

4 

0.0575 

4.90 

0.035 

8 

0.0446 

8.64 

0.07 

8 

0.0696 

9.26 

We  also  calculated  the  optical  properties  for  these  phantoms  (total  of  nine)  using  a  semi-infinite  fit  for  the  data  from 
a  strip  of  width  «  1.5  cm  at  the  top  part  of  the  scanning  region.  This  geometry  is  commonly  employed  and  was 
used  earlier  by  us  [4]  to  obtain  estimates  of  the  breast  optical  properties.  We  believe  that  the  main  drawback  of  the 
semi-infinite  model  is  that  it  ignores  the  multi-layer  boundary  structure  such  as  the  tissue  -Intralipid  boundaries.  It 
tends  to  underestimate  //s  and  both  /xa  and  y!s  “saturate”  at  higher  values  of  absorption  and  scattering. 

The  results  from  both  methods  are  shown  in  fig. (2)  where  we  plot  the  mean  calculated  properties  vs  the  expected 
values  from  multiple  measurements  for  each  phantom.  The  error  bars  are  derived  from  the  standard  deviations  from 
the  expected  values.  We  see  that  the  correlation  with  the  expected  values  for  the  finite  difference  method  is  within 
10%  whereas  for  the  semi-infinite  model  it  is  worse  than  45%. 

We  then  apply  this  method  for  the  data  from  healthy  volunteers.  Fig.  (3)  shows  the  histograms  of  optical  properties 
and  also  the  blood  saturations  for  seven  of  these  patients.  The  optical  properties  are  well  within  the  generally  accepted 
values.  The  blood  saturation  values  are  also  generally  satisfactory [6].  The  very  low  saturation  value  in  one  case  is 
probably  due  to  a  systematic  error. 

We  are  currently  acquiring  more  data  with  healthy  patients  and  we  still  have  a  large  data  set  obtained  that  needs  to 
be  analyzed.  The  model  is  also  being  tested  thoroughly  with  more  robust  phantoms.  This  work  will  be  presented  in 
April. 
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Fig.  1.  Models  used  for  finite  difference  calculations;  background  (left),  breast  (middle)  and  balloon  phantom 
(right)  are  shown  with  the  corresponding  boundaries.  The  shaded  regions  show  the  estimated  breast  volume 
in  the  segmentation  process. 


Fig.  2.  The  mean  calculated  absorption  (scattering)  coefficient  vs  the  expected  absorption  (scattering)  coef¬ 
ficient  is  shown  on  the  left  (right).  Solid  (dashed)  lines  show  the  finite  difference  (semi-infinite)  calculations. 
The  equality  line  (dotted)  is  also  shown.  The  error  bars  are  obtained  from  the  standard  deviations  from  the 
equality  line. 


Fig.  3.  Histograms  of  calculated  bulk  absorption,  scattering  coefficients  and  the  blood  saturation  of  the 
normal  volunteers’  breasts. 
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Abstract:  A  diffuse  optical  tomographic  imager  is  described.  It  employs  a  large  area  CCD  to 
collect  large  data  sets  (105  measurements)  from  32  sources  rapidly  (30s).  3D  images  of  breast 
tissue  phantoms  exhibit  -5mm  resolution. 
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1.  Introduction 

Diffuse  optical  tomographic  (DOT)  techniques  have  paved  the  way  for  several  groups  now  focusing  on  clinical 
research[l].  However  several  challenges  remain  for  developing  a  clinically  useful  tool.  One  of  the  main 
experimental  challenges  is  to  increase  the  speed  and  size  of  the  data  set  collected.  Coverage  and  resolution  demand 
larger  data  sets  and  yet  patient  movement  and  physiologic  noise  necessitate  fast  collection  speeds.  A  generally  held 
view  is  that  DOT  techniques  may  be  limited  to  between  2-10  mm  resolution  in  deep  breast  tissue.  One  might 
naively  assume  that  a  similar  sampling  density  would  be  sufficient.  For  typical  breast  volumes  (~1000cc)  this  would 
necessitate  between  103  (1cm  resolution)  and  105  (2mm  resolution)  measurements.  However  some  surface  features 
have  sub  millimeter  structure  and  are  near  the  high  sensitivity  regions  of  the  reconstruction  volumes.  For  these 
surface  structures  it  may  be  advantageous  data  set  sizes  of  greater  than  106.  Most  clinical  prototypes  thus  far  have 
much  smaller  data  sets,  i.e.  in  the  range  of  101  tolO3.  We  therefore  anticipate  the  potential  for  machines  with  larger 
data  sets  and  rapid  aquistion  to  substantially  improve  upon  existing  optical  mammography  devices. 

A  natural  tool  for  spatially  dense  light  sampling  is  a  CCD.  There  has  been  some  outstanding  work  in  the  biomedical 
optics  field  using  and  developing  CCD  systems  for  diffuse  light  diagnostics,  though  these  have  not  yet  applied  to 
clinical  mammography.  Time-resolved  (80psec  resolution)/2D-space-resolved  measurements  have  been  carried  out 
by  Alfano  and  coworkers[2,3].  Using  this  instrumentation  they  have  localized  a  fluorescent  object  in  chicken  breast, 
and  they  have  begun  to  work  on  the  general  image  reconstruction  problem  in  phantoms.  In  the  frequency  domain 
Sevick-Muraca's  group  has  adapted  homodyned  multi-pixel  RF  phase-sensitive  camera  technique  for  deep-tissue 
imaging  [4,5],  They  have  measured  images  of  DC/ AC  amplitude  and  phase  of  a  tissue  phantom  with  absorbing  and 
fluorescent  heterogeneities.  They  have  even  used  the  system  for  investigation  of  in-vivo  and  ex-vivo  animal  models 
with  mammary  cancer.  Boas  and  coworkers  have  also  recently  applied  CCDs  in  remission  geometry’s;  in  particular 
they  use  a  single  source  and  simply  measure  the  diffuse  light  amplitude  in  remission[6].  In  total  this  work,  and 
related  reconstruction  work  with  DC  data  [7],  demonstrates  the  feasibility  of  using  CCD’s  to  collect  data  and  as  well 
as  the  surprising  utility  of  DC  amplitude  only  data  for  image  reconstruction.  To  this  end  we  have  developed  clinical 
prototype  DOT  imager  using  a  larger  area  CCD  camera.  This  CW  imager  will  be  integrated  into  a  existing  slower 
scanning  RF  imager  that  we  are  currently  employing  to  study  bulk  breast  tissue  optical  properties. 

The  system  we  have  built  introduces  several  new  ideas.  Most  important  are  the  32  different  transmitted  image 
projections  of  the  illuminated  volume.  This  is  crucial  for  tomographic  analysis  of  the  data.  In  Figure  (1)  we 
illustrate  our  new  compression  plate,  intralipid  tank  and  CCD  imaging  system.  The  compression  plate  and  Intralipid 
tank  have  been  mechanically  redesigned  compared  to  our  previous  system.  In  particular  the  positional  repeatability 
and  to  stability  of  the  moving  compression  plate  relative  to  the  detection  window  were  improved.  In  addition  a 
modular  viewing  window  permits  the  use  of  diffuse  or  antireflective  transparent  windows.  Our  CCD  consists  of  a 
1300x1340  pixel  array  in  an  area  of  ~2.6cm  x  2.6  cm.  The  CCD  chip  is  read  by  a  16  bit  A/D  and  thermoelectrically 
cooled  (-40°C).  We  have  checked  the  light  transmission  through  a  range  of  tissue  phantoms  of  different  thicknesses 
in  order  to  ascertain  the  signal-to-noise  of  the  CCD  system.  We  chose  a  fairly  large  source-detector  plate  separation 
range  (~5-8cm),  and  made  an  Intralipid  tissue  phantom  with  (i,  0.05  cm-1  and  ps  =10  cm'1.  The  raw  CCD 
measurements  were  binned  in  8x8  square  units  to  give  an  equivalent  pixel  size  of  2mm  on  the  breast  surface.  Using 
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■  Figure  1.  CCD  imaging 

system  and  Intralipid  tank  with 
movable  compression  plate.  The 
modular  viewing  window  can  be 
either  a  diffuse  resin  window 
(p*=0.02crn !,  p^lOcm"1)  or  a 
transparent  window.  A  10-mW, 

809nm,  laser  light  source  is 
optically  switched  between  32 
fibers  arranged  in  a  hexagonal 
lattice  within  the  compression 
plate.  An  average  of  4mW  is 
delivered  to  the  fiber  tips.  A 
50mm  F#=l  .4  lens  images  the 
detection  plane  on  to  the  CCD 
chip.  Ambient  room  light  is 
eliminated  using  a  RG-9  optical 
glass  filter  and  light  shielding 
materials. 

known  detector  efficiency  and  gain  the  image  was  converted  to  photoelectron  values,  and  the  signal-to-noise  was 
assumed  to  scale  as  l/sqrt(Nphotoeiectrons).  From  our  clinical  measurements  we  know  that  most  breast  compressions  are 
~  6cm.  However  we  find  that  even  for  a  larger  breast  compression  of  7  cm  we  can  collect  ~3xl03  light 


measurements  with  a  S/N  of  better  than 
measurements  in  -30  seconds. 


A%  over  100cm  area.  With  32  available  sources  we  reach  a  total  of  10 


Figure  (2)  shows  the  results  from  a  measurement  made  with  two  (TA  and  TB)  absorbing  5x5x10  mm  targets  in  a 
suspension  of  Intralipid  with  a  source  detector  plane  distance  of  5.5cm.  The  source  plane  is  at  y=0  and  the  detector 
plane  is  at  y=  5.5cm  with  x-axis  oriented  horizontal  and  z-axis  oriented  vertical.  The  bulk  intralipid  tissue  phantom 
had  jia  =0.05  cm'1  and  jV  =10  cm'1  and  the  targets  (TA  ,TB)  had  (pa  =0.3  cm"1  ,pa  =0.2  cm"1)  respectively.  The  centers 
of  the  two  objects  were  located  in  the  middle  of  the  intralipid  tank  with  (Ta:x=-0.5,y=2.5,z=0)  and 
(Tb:x=1.5,y=2.5,z=0).  Ta  aligned  along  the  z-axis  and  Tb  aligned  along  the  x-axis.  To  demonstrate  the  tomographic 
capabilities  of  the  data,  3d  images  were  reconstructed  using  a  rytov  perturbation  approach  with  ART  matrix 
inversion.  For  these  preliminary  studies  computing  restrictions  necessitated  that  we  use  a  subset  of  the  available 
data.  In  particular  the  central  six  sources  were  used  located  at  z=0  and  spaced  evenly  at  1.4cm  steps  in  the  x 
direction.  The  detectors  were  then  binned  to  4. 6x4.6  cm  area  measurements  covering  a  -4x1 0cm  field  of  view  with 
19  positions  in  x  and  1 1  positions  in  z.  A  total  of  1254  measurements  were  available.  From  this  total  detector  set  a 
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Figure  2.  The  central 
6x2. 5x2. 5  volume  (x,y,z)of 
the  entire  10x5.5x4cm 
volume.  Two  slices  of  the 
3D  pa  image  maps  are 
shown,  a)xz  plane  at 
y=2.5cm  for  the  phantom 
b)xy  plane  at  z=0  for  the 
phantom,  c)xz  plane  at 
y=2.5cm  for  the 
reconstrution  d)xy  plane  at 
z=0  for  the  reconstrution,  e) 
e)Cross  sections  at  y=0,z=0. 
F)  cross  section  at  x=- 
0.5, z=0. 


reduced  data  set  was  selected  using  a  noise  threshold  defined  relative  to  the  peak  light  intensity  for  each  source. 
Ultimately  1001  total  measurements  were  used  for  the  reconstruction’s.  Using  rytov  absorption  weights,  the 
10x5.5x4cm  (x,y,z)volume  was  meshed  in  2.3x2.3x2.3mm  voxels.  An  Algebraic  Reconstruction  Technique  (ART) 
was  used  to  invert  the  measurements  an  obtain  the  3D  pa  images.  Figure  2  shows  the  central  slices  in  the  xz  plane  (at 
y=0)  and  the  xy  plane  (at  z=0)  for  the  (a&b)  actual  pa  values  and  (c&d)  the  reconstructed  values.  For  quantitative 
comparison  the  1  dimensional  cross  sections  e)  in  the  x  direction  at  y=0,z=0  and  f)in  the  y  direction  at  x=-0.5,z=0 
are  shown.  Note  that  the  orientations  of  the  two  objects  in  the  xz  plane  are  observable.  Thus  the  resolution  in  the  xz 
plane  is  ~5mm.  The  reconstructed  values  are  ~2  times  smaller  than  expected.  This  can  be  attributed  to  the  reduced 
(~lcm)  resolution  in  the  y  direction  that  spreads  the  optical  weight  over  a  larger  volume  leaving  the  total  integrated 
signal  8j.ia*volume  the  same.  It  is  anticipated  that  reconstruction’s  that  utilize  more  of  the  available  source  views 
will  improve  on  the  resolution  in  the  y  direction  and  correspondingly  the  quantification  of  the  optical  properties. 
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Figure  3.  Image  error  (r)  as  a  function  of  the  number  of 
measurements  used  in  the3D  reconstruction. 
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To  test  the  effect  of  data  size  an  image  quality 
criteria  is  needed.  For  our  comparison  we  use  a 
standard  image  error  function  (ref)  applied  to 
the  central  xz  plane.  In  particular  we 
define 

Y  X  Y  X 

r=  Z  Z  1 1.  ,-r  |/  Z  I  I'- 

/  =  1  y  =  1  ,J  i  =  lj  =  1 

where  X,  Y  are  the  dimension  in  voxels  of  the 
reconstructed  image,  t  is  the  expected  optical 
property,  and  r  is  the  reconstructed  optical 
property.  The  data  size  is  altered  by  sparsely 
selecting  the  pixels  while  leaving  the  field  of 
view  constant.  Figure  3  shows  the  error  as 
function  of  data  set  size.  Here  the  true  noise 
threshold  truncated  data  size  is  used.  Note  that 


we  are  continuing  to  improve  the  image  quality  as  we  increase  are  detection  density  to  4.6  mm.  It  is  expected  that 
with  more  heterogeneous  samples  the  effects  will  be  even  more  pronounced.  Future  work  will  extend  this  analysis  to 
larger  data  sets  and  more  complex  tissue  phantoms. 


In  conclusion  we  have  developed  a  clinical  prototype  DOT  imager  using  a  larger  area  CCD  camera  that  records  105 
measurements  in  30  seconds.  This  is  a  significantly  larger  data  set  obtained  more  quickly  than  typical  existing 
clinical  optical  mammography  machines.  Preliminary  3D  reconstructions  of  tissue  phantom  demonstrate  -5mm 
resolution  with  a  limited  data  set.  We  will  integrate  this  system  into  a  existing  slower  scanning  RF  imager  that  we 
are  currently  employing  in  the  clinic. 
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Abstract:  A  new  integro-differential  equation  for  diffuse  photon  density 
waves  (DPDW)  is  derived  within  the  diffusion  approximation.  The  new 
equation  applies  to  inhomogeneous  bounded  turbid  media.  Interestingly,  it 
does  not  contain  any  terms  involving  gradients  of  the  light  diffusion 
coefficient.  The  integro-differential  equation  for  diffusive  waves  is  used  to 
develop  a  3D-slice  imaging  algorithm  based  on  the  angular  spectrum 
representation  in  the  parallel  plate  geometry.  The  algorithm  may  be  useful 
for  near  infrared  optical  imaging  of  breast  tissue,  and  is  applicable  to  other 
diagnostics  such  as  ultrasound  and  microwave  imaging. 
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1.  Introduction 

Tomographic  imaging  of  deep  tissues  with  diffusive  waves  makes  possible  functional 
imaging  based  on  novel  optical  contrasts  derived  from  tissue  chromophores,  structure  and 
metabolism  [1-3].  The  near-infrared  spectroscopic  properties  of  breast  tumors  for  example, 
have  been  found  to  differ  from  adjacent  normal  tissues.  Such  spectroscopic  signatures  hold 
promise  for  increased  tumor  detection  sensitivity  and  specificity  [4-10].  In  order  to  take  full 
advantage  of  these  new  contrasts  it  is  critical  to  develop  high  fidelity  three-dimensional 
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optical  imaging  algorithms  for  diffusing  light  in  turbid  media  such  as  the  breast.  To  this  end  a 
number  of  image  reconstruction  schemes  have  been  developed  [11-13].  Many  of  these 
schemes  use  integral  equations  based  on  Born  or  Rytov  approximations  to  generate  a  set  of 
linear  equations  which  are  then  solved  to  update  the  absorption  and  scattering  coefficients 
associated  with  each  voxel  in  the  reconstructed  volume.  The  integral  formulation  is  attractive 
because  of  its  speed.  However,  to  our  knowledge  most  of  these  formulations  ignore  terms 
involving  gradients  of  the  light  diffusion  coefficient  [14,15]  Although  this  approximation  is 
often  reasonably  accurate,  a  recent  paper  [17]  has  established  the  importance  of  this  gradient 
term,  showing  that  its  neglect  is  responsible  for  cross  coupling  between  absorption  and 
scattering  images. 

In  this  paper,  we  derive  a  new  integro-differential  equation  for  diffuse  photon  density 
waves  (DPDW)  within  the  diffusion  approximation.  The  new  equation,  which  is  developed 
for  bounded  turbid  media  does  not  explicitly  contain  any  terms  involving  gradients  of  the 
light  diffusion  coefficient.  We  then  use  this  integro-differential  equation  and  develop  a 
theoretical  inverse  scattering  algorithm  for  three-dimensional  image  reconstruction.  The 
theoretical  framework  follows  the  principles  of  near-field  diffraction  tomography  based  on 
the  angular-spectrum  representation;  in  particular  we  show  how  to  develop  three-dimensional 
slice  images  of  a  breast  compressed  between  two  parallel  plates.  The  technique  employs  a 
series  of  plane  diffuse  photon  density  waves  with  different  modulation  frequencies.  For  this 
set  of  plane  incident  diffusive  waves,  the  algorithm  requires  two-dimensional  FFT  operations, 
and  a  one-dimensional  matrix  inversion.  Thus  the  method  is  non-iterative  in  two  space 
dimensions  and  is  therefore  computationally  fast. 

2.  Integro-differential  equation  for  diffuse  photon  density  waves 

We  begin  with  the  equations  satisfied  by  the  diffuse  photon  density,  O  and  the  diffuse 
photon  flux  density,  J  [18] 


(1) 

3 £(?•<») 

(2) 

(--+ m.M 

c 

Here,  0)  denotes  the  modulation  frequency  of  the  source  intensity,  // a  and  jU  *s  are  the 
absorption  and  transport  scattering  coefficients  of  the  bounded  turbid  media,  c  is  the  velocity 
of  light  in  the  medium  and  S0  (r ,  0)) ,  S{  (r ,  CO)  are  the  monopole  and  dipole  moments  of  the 

intensity  modulated  optical  source  respectively.  From  these  equations  we  eliminate  the 
photon  flux  density  and  obtain,  after  a  lengthy  calculation,  the  following  wave  equation 
involving  only  the  photon  density  function 
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Here 


kd2  =3jU's(--Ma)  (4) 

c 

is  the  square  of  the  complex  wavenumber  kd .  The  bar  over  a  physical  parameter  denotes 

spatial  average  of  that  parameter.  Notice  that  for  unmodulated  source  intensity  (the 
continuous  dc  domain  in  contrast  to  the  frequency  domain),  the  photon  density  is  a  non¬ 
propagating  damped  wave.  The  wave  number  in  this  case  is  purely  imaginary.  For  frequency 
modulated  source  intensity,  the  wave  number  is  complex  and  the  photon  density  acquires  the 
characteristics  of  a  damped  or  attenuated  wave. 

Because  of  its  differential  nature  Eq.  (3)  will  admit  many  solutions  and  boundary 
conditions  are  necessary  to  obtain  the  physical  solution.  The  problem  of  finding  the  proper 
boundary  condition  is  complex  and  is  generally  further  exacerbated  because  the  diffusion 
approximation  is  not  valid  near  the  boundary.  Inspite  of  this  basic  difficulty,  diffusion  theory 
is  used  widely  due  to  its  simplicity  and  the  applicability  of  numerous  analytical  and  numerical 
techniques  for  its  solution  in  finite  domains  of  arbitrary  geometry.  Approximate  boundary 
conditions  are  generally  adopted  which  have  been  shown  to  be  fairly  accurate.  Here,  we  adopt 
the  general  boundary  condition  [19] 


h-J(r,C0)  =  0&(r,C0). 


(5) 


The  value  of  a  can  be  defined  based  on  boundary  considerations  or  taken  as  a  fitting 
parameter  for  a  given  interface.  In  Eq.  (5)  h  is  the  unit  outward  normal  to  the  boundary 
surface.  Using  Eq.  (5),  it  is  now  possible  to  transform  Eq.  (3)  into  an  integral  equation,  the 
solution  of  which  will  provide  one  with  a  physical  picture  of  the  interaction  of  the  diffuse 
photon  density  wave  with  the  turbid  media.  The  basic  interaction  of  light  with  the  molecules 
of  the  turbid  media  is  already  included  in  the  absorption  and  scattering  parameters.  The 
determination  of  these  material  parameters  is  at  the  heart  of  the  optical  modality  for  cancer 
diagnostics.  The  integral  representation  forms  a  basis  for  extraction  of  these  parameters  from 
measurements  of  the  diffusion  waves  at  the  boundary.  We  next  derive  the  integro-differential 
equation  using  Eqs.  (3)  and  (5). 

We  follow  the  standard  Green's  function  technique  for  the  derivation  of  an  integral 
equation  from  a  partial  differential  wave  equation  [20].  However,  we  deliberately  choose  a 
simple  Green’s  function  corresponding  to  the  operator  on  the  left  hand  side  (L.H.S.)  of  Eq. 
(3)  that  satisfies  the  equation 

(V2  +  k/)G(r,  r\  CO)  =  -4 7lS(j  -  r\  (6) 


and  is  given  by 
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with 


(8) 


Here  lm(kd  )  >  0.  Notice  that  our  choice  of  an  infinite  space  Green's  function  will  not  limit 
in  any  way  when  we  consider  finite  domains;  finite  medium  effects  are  properly  treated  by 
the  boundary  integrals.  If  we  now  multiply  both  sides  of  Eq.  (3)  and  (6)  by  Gd  (r ,  r  ,  OS)  and 

(J)(r,  CO)  respectively  and  subtract  the  resulting  equations  from  one  another  and  use  a 
number  of  standard  Green’s  identities,  we  obtain 


f S0(r\co)Gd(r,r\a)d3r'  +  -^  [ Sl(r\a).V'Gd(r,r\(0)d'*r' 
4  7T  Jv  Aft  * 

+  —  J(3  ajusGd(r,  r',0) )  +  V'Gd(r,r',  a>)d,S' 

4  ft  s 

4  7t 


+  —  f  a>)YGd  (r\ r)d 3r. 

4xJv  jus{7') 


(9) 


In  deriving  Eq.  (9),  we  also  made  use  of  Eqs.  (l)-(2)  and  (5).  V  denotes  the  volume  of  the 
turbid  medium  and  the  surface  bounding  the  volume  is  denoted  by  S.  The  first  two  terms  on 
the  right  hand  side  (R.H.S.)  of  Eq.(9)  are  related  directly  to  the  source  and  the  second  term  on 
the  R.H.S.  of  Eq.(9)  arises  because  of  the  presence  of  the  boundary.  The  solution  to  this 
integro-differential  equation  is  unique  as  long  as  there  exists  no  non-trivial  solutions  to  the 
corresponding  homogeneous  integral  equation.  It  is  possible  that  for  certain  modulation 
frequencies  the  resulting  homogeneous  integral  equation  may  possess  non-trivial  solutions. 
These  frequencies  will  generally  be  complex  and  will  correspond  to  resonance  behavior  of  the 
photon  density  waves.  We  shall  not  dwell  with  this  singular  case,  but  refer  the  reader  to  a 
related  paper  on  resonance  theory  by  one  of  the  authors  [15].  Our  Eq.(9)  differs  from  the 
corresponding  equation  used  in  reference  [14],  in  which  the  boundary  effect  was  ignored.  The 
treatment  of  the  source  contribution  to  the  integral  equation  is  also  different. 

For  a  homogeneous  diffusing  media,  the  integro-differential  equation  becomes 
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Equations  (9)  and  (10)  are  the  basic  equations  of  our  paper.  They  are  derived  using  the 
diffusion  approximation  to  photon  transport  and  without  any  other  approximations.  Eq.  (9)  is 
an  intro-differential  equation  for  the  photon  density  wave  and  involves  only  the  absorption 
and  transport  scattering  coefficients.  Notice  that  Eq.  (9)  does  not  explicitly  contain  gradients 
of  light  diffusion  coefficient.  Our  derivation  is  exact  and  does  not  assume  any  condition  on 
the  gradient  of  the  light  diffusion  coefficient.  As  has  been  pointed  out  recently  [17]  neglect  of 
the  gradient  term  leads  to  cross  coupling  between  scattering  and  absorption  images.  Eq.  (10) 
describes  the  homogeneous  or  the  background  photon  density  field  that  includes  the  source 
contributions  (the  first  two  terms  on  the  R.H.S.),  as  well  as  contributions  from  the  boundary 
(the  last  term  of  the  R.H.S.).  Note  that  Eq.  (10)  contains  only  average  values  of  the  scattering 
parameters  and  the  surface  integral  includes  the  total  field. 

Eq.  (10)  relates  the  value  of  the  diffuse  photon  density  wave  at  any  point  inside  the  turbid 
medium  to  the  monopole  and  dipole  moment  of  the  source  and  to  the  photon  density  on  the 
surface  bounding  the  medium.  A  surface  integral  equation  for  the  boundary  value  of  the 
photon  density  is  obtained  if  we  take  the  point  7  to  be  on  the  surface.  This  limiting  operation 
must  be  done  carefully  because  the  integrals  involving  the  gradient  of  the  Green's  function  are 
not  continuous  when  the  interior  point  approaches  a  point  on  the  surface.  The  other  integrals 
involving  the  Green's  function  however  are  continuous  [21].  From  an  experimental  point  of 
view,  if  the  photon  density  on  the  surface  surrounding  the  turbid  medium  is  measured,  then 
one  can  use  Eq.  (10)  to  determine  the  diffuse  photon  density  wave  anywhere  inside  the 
medium  by  simply  carrying  out  the  surface  integrals  in  Eq.  (10).  If  the  value  of  this  estimated 
photon  density  matches  that  of  the  measured  value  at  all  points  on  the  surface,  then  an 
absence  of  tissue  heterogeneity  is  indicated.  In  practice  there  will  never  be  perfect 
cancellation,  but  a  larger  difference  indicates  the  presence  of  stronger  heterogeneities  inside 
the  tissue. 

Another  approximate  way  to  evaluate  the  background  field  is  to  evaluate  the  surface 
integrals  in  Eq.  (10)  iteratively.  First  ignore  the  third  surface  term  involving  the  boundary 
values  of  the  photon  density  in  Eq.  (10)  and  use  only  the  first  two  terms  involving  the  source 
terms  to  calculate  the  first  order  background  field  everywhere.  Then  update  the  value  of  the 
background  field  by  using  the  first  order  values  in  the  last  integral  in  Eq.  (10).  It  is  important 
to  devise  a  way  to  estimate  a  reasonable  background  field,  because  when  this  field  is 
subtracted  from  the  total  measured  diffuse  photon  density  field,  then  the  remainder  is  the  field 
directly  attributable  to  the  presence  of  the  heterogeneity. 

We  have  already  noted  that  unlike  other  integral  formulations  our  integro-differential 
equation  does  not  contain  an  explicit  term  involving  the  gradient  of  the  transport  scattering 
coefficient.  This  is  important  in  the  inverse  problem  where  the  absorption  and  transport 
scattering  heterogeneity  are  to  be  estimated  from  the  measured  values  of  the  field.  The 
formalisms  that  include  the  gradient  term  will  of  necessity  have  to  use  a  numerical  estimate  of 
the  gradient  value.  Such  procedures  may  lead  to  numerical  instabilities  and  numerical 
artifacts.  We  next  develop  a  3D-slice  image  reconstruction  technique  based  on  the  integro- 
differential  equation  for  the  photon  density,  i.e.  theEqs.  (9)  &  (10). 


3.Inversion  algorithm  for  3D-slice  imaging 


We  make  the  assumption  that  the  heterogeneities  are  weak  and  do  not  perturb  the  background 
field  substantially.  Then  we  can  approximate  the  photon  density  within  the  integral  operator 
with  the  background  field.  With  this  first  order  Born  approximation,  Eq.  (9)  becomes 

Op(r,  co)  -  (7, co)  - o"  (7, CO)  =  ( 7')-Jia )0>w  (?',  co)Gd(7\  7, co)d 3 

t\7l  y 


1 

+  - 

4  71 


f  )VV"  (7,  co)YGd  (7, 7,  co)d V. 

v  Ms(r  ) 


(11) 
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The  superscript  M  implies  measured  field  and  the  superscript  P  denotes  the  processed 
(scattered)  field  values.  In  Eq.(12),  we  have  used  the  identity 


<D 


3  X 


H (F,  CD)  =  I*  S0(r',  co)Gd (F, r', co)dir'  +  -—  f  5, (F', co)AJ'Gd (Jr,  r\ co)d3r' 
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+  —  J  (3  ajtsGd  (r,  r',  CD)  +  V'Gd  (F,  r',  co).h)<t>M  (r\  co)dS'.  (12) 

Aft  s 

The  quantity  (F,  CO)  often  identified  as  the  background  field  can  be  determined  from  the 
knowledge  of  the  source  and  the  measured  field  at  the  boundary  from  Eq.  (12). 

Although  our  integro-differential  equation  is  valid  for  arbitrary  geometry,  from  now  on 
we  shall  restrict  ourselves  to  the  case  of  a  plane  parallel  plate  configuration.  The  plane 
parallel  plate  is  an  accepted  clinical  diagnostic  configuration  for  breast  imaging.  It  offers  the 
possibility  of  compressing  the  breast  to  decrease  its  thickness  for  increased  signal  strength. 
We  also  take  the  transverse  dimension  of  the  plate  to  be  sufficiently  larger  than  the  plate 
separation.  Because  of  the  negligibly  small  value  of  the  field  at  the  periphery  of  sufficiently 
large  plates,  in  practice  the  transverse  plate  extent  can  be  taken  to  be  infinite.  In  essence,  Eq. 

(11)  is  an  integral  equation  for  jUa(r)  ,  and  //'(F)  with  all  the  other  variables  such  as  the 
background  field,  the  Green's  function  and  the  background  optical  parameters  either  given  or 
estimated. 

In  the  new  integro-differential  equation,  we  use  the  Weyl  representation  [22-23]  of 
Green's  function 


G(r,r',00)-  \dp  (*  dq—^—exp(ip(x-  x')  +  iq(y  -y')  +  im\z~  z\), 
J  J  2ftm 


(13) 


where  m  =  <Jkd2  -  p2  -  q2  ,  with  Im(m)  >  0. 

The  direction  z  is  chosen  perpendicular  to  the  surface  normal  and  the  x  and  y  directions  are 
chosen  along  the  transverse  directions.  There  is  a  singularity  in  the  z  direction  as  can  be  seen 
from  the  absolute  value  operator  in  the  exponent  of  the  Weyl  representation.  This 
representation  is  also  known  as  the  angular  spectrum  representation  and  is  fundamental  to  the 
theory  and  algorithm  now  commonly  known  as  Diffraction  tomography  [24-25].  Diffraction 
tomography  was  applied  to  far  field  coherent  optical  imaging  systems,  but  recently  this 
technique  has  been  extended  to  the  near  field  diffuse  photon  density  wave  based  imaging  [26- 
34].  With  this  representation  of  Green's  function  the  integrals  in  x  and  y  become  convolution 
integrals  which  turn  into  algebraic  products  in  p,  q  space.  Taking  the  two-dimensional 

Fourier  transform  of  Eq.  (12)  at  the  detector  plane  z  =  Zd  we  find  that 


®\p,q  ,Zd,CO)  =  ]dz'Ta(p,q,z  ,CO)K(p,q,zd  ~  Z  ,CD) 

Zs 

Zd  -  , 

+  J  dz'Ts  (p,  q,  z',  CD)  ■  kK(p,  q,  zd  ~z,co), 

Zs 


(14) 
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where 


and 


Ta(r ,ct))  -  (/4(0-/4)O  ( r,(o ), 

4  x 

(15) 

f,  (r,  CO)  =  -F  (1  -  V  0H  (F,  CO), 

(16) 

4sr  D 

k{p,q,zd  z  , of)  —  exp (im(zd~z)). 

2mn 

(17) 

In  Eqs.(14)-(17),  k  —  xp  +  yq  +  zm  is  the  wave  vector  and  Zs,Zd  denote  the  positions  of 

the  source  and  detector  plates  respectively.  The  caret  (  A)  over  functions  signifies  the  Fourier 
transform  of  a  function.  Several  things  are  worth  noting  at  this  point.  First,  the  absorption 
coefficient  contribution  depends  on  a  scalar  function,  while  the  transport  scattering  coefficient 
depends  on  the  product  of  two  vectors.  Second,  the  gradient  of  the  diffuse  photon  density 
wave  interacts  with  the  scattering  heterogeneity,  while  the  diffuse  photon  density  wave 
interacts  with  the  absorption  heterogeneity.  Finally,  the  inverse  of  the  transport  scattering 


coefficient  has  been  explicitly  replaced  by  the  photon  diffusion  coefficient  D  — 


c 


3a; 


Eq.  (14)  is  a  very  simple  looking  equation,  but  in  fact  a  lot  of  physics  is  hidden.  The 
Fourier  transform  of  the  so-called  tumor  functions,  Ta  and  Ts  which  are  products  of  the 

scattering  parameters  and  the  background  field  in  real  space,  will  be  convolution  integrals  in 
the  Fourier  space.  This  convolution  in  the  spatial  frequency  space  between  the  tumor 
parameters  and  the  background  field  couples  temporal  frequency,  CO  with  spatial  frequency 
(p,q).  This  implies  that  the  use  of  a  set  of  modulation  frequencies  does  not  in  any  easy  manner 
let  us  solve  for  the  spatial  Fourier  amplitudes  of  the  tumor  functions.  This  coupling  is 
eliminated  if  the  background  field  is  spatially  in  the  form  of  a  plane  wave  rather  than  a 
spherical  wave.  The  usefulness  of  plane  wave  sources  is  well  known  in  the  field  of  diffraction 
tomography  [24]. 

We  will  use  a  series  of  plane  waves  whose  central  modulation  frequencies  are  stepped  to 
span  some  frequency  bandwidth.  For  the  plane  wave  case  the  tumor  functions,  fa  ,  Ts ,  can 


be  written  as  a  product  of  two  functions.  One  function,  ta  or  tD  ,  contains  information  about 

the  heterogeneities  and  depends  upon  (p,  q,  z)  and  the  other  function  depends  upon  the 
background  plane  waves  and  modulation  frequency  (see  Eq.  (18)  below).  One  can  then  utilize 
this  temporal  frequency  bandwidth  to  obtain  resolution  in  the  longitudinal  (z)  direction.  The 
key  advantage  of  using  the  series  of  stepped  central  modulation  frequencies  is  that  the 
detection  system  is  still  narrow  band  and  one  can  use  maximum  permissible  signal  strength  at 
all  discrete  frequencies.  We  plan  to  publish  in  a  later  paper  the  details  of  this  technique  with 
simulated  data,  but  here  we  point  out  the  salient  features  of  the  method.  For  the  plane  wave 
case  we  find 


Tu(p,q,z,(D)  =  ta(p,qtz)[®  +  exp(  +im0z)  +  O  "  exp(  - imQz )}, 
fs{p,q,z,a >)  =  td(p,q,z){ka+®  +  exp(  +im0z)  +  k0~<&~  exp(  -imQz)}, 
ia(P’<hz)  =  -j-\.Mad(P  ~  Pn)d(q  ~  <la )~  Mu(P  ~  Po.<l  ~  <?<>)]> 

4  71 


#9231  -  $15.00  US 
(C)  1999  OS  A 


Received  March  05,  1999;  Revised  April  09,  1999 
12  April  1999  /  Vol.  4,  No.  8  /  OPTICS  EXPRESS  238 


and 


(18) 


tD (p,  q,  z)  =  i  [D  S{p  -  p0  )S(q  -  q0 )  -  D(p  -  p0  q  -  q0 )] , 
k0+  =  Wo  +  Wo  +  zm0 ,  k~  =  xp0  +  yq0  -  zm0. 


Here  p0  -  Re(kd)cos(0)cos(0)  ,q0  =  Re(kd)cos(0)sin(0)  where  0  is  the  angle 

the  plane  wave  makes  with  respect  to  the  z  direction,  and  (j)  is  the  corresponding  azimuthal 
angle.  Care  must  be  exercised  in  the  angle  interpretations  as  we  are  dealing  with  attenuating 
plane  waves.  For  the  case  of  a  normal  incidence,  p0=  0  and  q0  =  0 .  The  parameter 

m0  =  Jk/  -  Po  ~  Po  with  Im(m0)  >  O.The  general  mode  structure  of  the 

background  plane  wave  for  the  parallel  plate  slab  geometry  can  be  shown  from  the  solution  of 
the  homogeneous  equation  ( Eq.  (10))  to  be  a  sum  of  both  forward  and  backward  propagating 
waves  [33] 

®0(x,  y,  z,  cd. )  =  ®+  exp (ip0x  +  iq0y  +  im0z )  +  O-  exp (ip0x  +  iq0y  -  im0z).  (19) 

Here,  e>+,  and  ®  are  complex  constants  obtained  from  solution  or  measurements  as 
discussed  earlier  in  the  section  on  the  background  field. 

We  approximate  the  integration  in  Eq.  (14)  in  the  z  variable  by  means  of  a  sum  over  N 

slices  between  the  parallel  plates  (the  slice  thickness  can  be  different  and  is  denoted  by  A Zj 
with  j  ranging  from  1 :  N),  and  use  Eqs.  (15)-(19)  to  obtain  the  following  result 

N 

®(zd,cOj)  =  ^J[ta(zi)Mzi,(Oj)  +  id(zi)f2(zi,0)j)]  for  j  =  l:2N,  (20) 
1=1 


where 


/i  (z,  >  coj )  =  A  Zj  {<E>+  exp(+;m0-, )  +  O"  exp(-/m0z,-  )}K(zd  -  z,- ,  co), 
fi  ( Zt  ,(Oj)  =  -Az,  { (ppQ  +  qq0  +  mm0)<Z>+  exp(+tm0z, )  + 

( pp0  +  mo  ~  exp(-/m0z, )  }K(zd  -  z, ,  co).  (21) 


Here,  we  suppressed  the  (p,q)  arguments  for  the  sake  of  brevity.  Eq.  (20)  provides  2N  linearly 
independent  complex  equations  that  can  be  solved  for  the  two  dimensional  Fourier  transform 
of  the  absorption  and  scattering  coefficients  at  the  N  selected  slice  positions  in  the  z  direction. 
This  can  be  seen  more  clearly  when  Eq.  (20)  is  written  in  a  matrix  form 

-  =  £‘-  (22) 


Where  Ois  the  vector  of  2N  dimension  consisting  of  the  measurements  of  the  processed 
(scattered)  field  at  the  2N  modulation  frequencies,  /  is  a  2NX2N  matrix  consisting  of  two 
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2NXN  matrices  /,  ,  and  f2  defined  in  Eqs.  (21a)-(21b).  Eq.  (22)  can  be  solved  by  direct 
matrix  inversion  as  shown  below  or  iteratively  by  following  methods  such  as  SIRT,  ART  [13] 
.  f  consists  of  the  sum  of  the  two  N  dimensional  vectors  ta  and  tD. 


i=Or'& 

l  in  Eq.  (23)  is  the  two-dimensional  Fourier  transform  of  the  absorption  and  scattering 
coefficients  at  each  of  the  N  number  of  slices  in  the  direction  normal  to  the  plates  (the  z 
direction).  The  values  of  the  scattering  coefficients  in  the  direction  transverse  to  the  surface 
normal  (the  x  and  y  directions)  are  obtained  by  applying  the  two-dimensional  inverse  Fourier 
transform,  i.e. 

^  oo  oo 

t(x, y,zi)  =  — —  f  f dpdq(f(p, q))~'  •  q) exp (-ipx - iqy )  for  i  =  l:N 

\7T  j  j  = 


For  the  Breast  imaging  case,  Eq.  (24)  provides  one  with  3D-slice  images  of  the  absorption 
and  scattering  coefficients  at  the  selected  "virtual"  slice  positions  in  the  interior  of  the  breast 
compressed  between  the  two  plates.  These  images  may  be  reconstructed  with  light  of 
different  wavelengths  within  the  near-  infrared  window.  Such  measurements  improve  the 
usefulness  of  nearfield  optical  imaging  modality  employing  diffuse  photon  density  waves. 

4.Conclusion 

We  presented  a  rigorous  derivation  of  a  new  integro-differential  equation  for  the  diffuse 
photon  density  wave  in  an  inhomogeneous  bounded  turbid  medium  within  the  diffusion 
approximation.  The  resulting  equation  contains  only  the  absorption  and  scattering  coefficients 
and  does  not  include  gradients  of  the  light  diffusion  coefficient.  We  also  derived  a  novel 
surface  integral  equation  for  the  diffuse  photon  density  waves  for  the  case  of  a  bounded 
homogeneous  turbid  medium.  We  used  the  theory  to  develop  a  detailed  imaging  algorithm 
that  provides  one  with  3D-slice  images  of  the  absorption  and  scattering  coefficients  of  tissue 
(breast)  that  is  compressed  between  two  parallel  plates.  We  believe  this  algorithm  will  be 
useful  in  increasing  the  specificity  of  the  near  infrared  optical  imaging  modality.  Our 
algorithm  being  based  on  inhomogeneous  wave  equation  of  the  Helmholtz  type  may  also  be 
applicable  to  other  imaging  modalities  such  as  ultrasound. 
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Abstract:  Three-dimensional  diffuse  optical  tomography  is  a  computationally  intensive 
procedure.  An  iterative  perturbative  algorithm  is  followed  and  parallelized.  The  results  show  that 
three-dimensional  images  can  be  reconstructed  on  realistic  time  scales  compared  to  the 
synchronous  approach. 
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OCIS  codes:  (100.6890)  Three-dimensional  image  processing;  (170.6960)  Tomography;  (170.5280)  Photon  migration 


1.  Introduction 

Most  Diffuse  Optical  Tomography  (DOT)  has  occurred  in  2-D.  In  the  best  of  these  cases,  a  cylindrical  geometry 
was  employed  to  reduce  the  dimensionality  of  the  problem  [1,2].  In  fact,  only  recently  have  full  3-D  reconstructions 
been  accomplished  and  compared  to  2-D  [3,4];  not  surprisingly,  the  3-D  images  were  superior.  Additionally,  in  any 
clinical  setting  the  problem  is  necessarily  3-D.  Therefore,  we  are  convinced  that  3-D  reconstruction  is  the  best 
approach  for  optimal  image  fidelity  in  a  clinical  setting. 

The  main  barriers  for  3-D  reconstruction  are  memory  requirements  and  computational  processing  time.  Arridge 
and  Schweiger  have  presented  an  algorithm  that  reduces  memory  requirements  significantly  [5].  Unfortunately,  the 
computational  processing  time  was  still  significant.  Subsequently,  they  recently  presented  a  3-d  segmented 
reconstruction  using  parallel  processing  [4].  In  this  application  they  parallelized  the  forward  problem  matrix  solver 
and  received  a  factor  of  four  speed-up  when  using  eight  processors.  While,  this  is  promising  it  is  still  necessary  to 
speed  up  the  reconstruction  time  as  much  as  possible  in  order  to  make  3-d  clinical  imaging  a  reality.  Therefore,  we 
have  developed  a  parallel  3-d  DOT  algorithm  that  parallelizes  the  algorithm  itself  and  plans  for  expanding  the 
parallelization  are  underway. 

2.  Methods 

A  variety  of  methods  have  been  developed  for  (DOT).  These  include  fitting  to  an  analytic  solution  [6-8], 
backprojection  method  [9-11],  diffraction  tomography  in  k-space  [12-16],  perturbation  approach  [17-23],  elliptic 
systems  method  (ESM)  [24,25],  and  direct  method  [26].  Each  method  has  advantages  and  drawbacks,  but  the  only 
method  to  date  that  has  been  utilized  in  full  3-D  simultaneous  reconstructions  of  absorption  and  scattering 
coefficients  for  complicated  geometries  and  boundary  conditions  is  the  perturbation  approach.  Within  the 
perturbation  approach  there  are  a  couple  of  methods  to  solving  the  forward  problem  and  determining  the  Jacobian 
[for  a  review  see  27].  Additionally,  there  are  a  variety  of  methods  for  solving  the  inverse  problem  [2,5,27-31]. 

Our  algorithm  is  based  on  an  iterative  perturbative  approach  that  uses  the  explicit  adjoint  formulation  for  the 
inverse  problem  and  parallel  processing  in  master-slave  style  (fig.  1).  This  requires  solving  the  forward  problem  for 
each  source  (NS=number  of  sources),  solving  the  adjoint  problem  for  each  detector  (ND=number  of  detectors),  and 
solving  the  inverse  problem  ([NMxNV],  NM=number  of  measurements,  NV=number  of  voxels).  In  the  first  step 
each  node  has  its  own  set  of  arrays  and  variables  that  are  initialized.  During  the  second  step  the  master  sends  out  the 
optical  properties,  which  in  the  first  iteration  are  just  the  background  values.  In  the  third  step,  the  slaves  calculate 
the  solutions  to  the  forward  problem;  consequently  NS  slaves  are  used  in  this  step.  The  solution  vectors  are  then 
returned  to  the  master  and  at  this  point  the  master  checks  for  convergence.  If  convergence  has  not  be  achieved  then 
the  slaves  calculate  the  solutions  to  the  Green’s  function  resulting  in  ND  slaves  being  utilized.  Again  the  solution 
vectors  are  returned  to  the  master  and  the  Jacobian  is  determined.  Next  the  master  solves  the  inverse  problem  by 
using  a  spatially  variant  regularized  conjugate  gradient  optimization  method  [2,21].  Finally,  the  optical  properties 
are  updated  on  the  master  and  the  algorithm  repeats  until  convergence  has  been  achieved. 

The  outlined  approach  utilizing  parallel  processing  results  in  many  gains.  First  we  are  able  to  calculate  our 
forward  solutions  simultaneously.  Therefore,  what  would  normally  take  NS  x  Atforward  (Atforward  =  the  time  of  one 
forward  solve)  is  now  {NS/Min(NS,  NP  (number  of  processors))}  x  Atforward  plus  some  small  amount  of  time  for 
communication.  Similarly,  for  the  Green’s  function,  processing  is  reduced  from  ND  x  AtGreen  (AtGreen  =  the  time  of 
one  Green’s  function  solve)  to  {ND/Min(ND,  NP  (number  of  processors))}  x  AtGreen-  If  NP>NS,ND  then  our  matrix 
solve  time  is  reduced  by  ~  NS+ND. 
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Fig.  1 :  Algorithm  for  full  3-D  reconstruction. 
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Fig.  2:  Maps  of  A|Lia(left)  and  A|V  (right)  actual  and  reconstructed. 


3.  Results 

A  full  3-D  parallel  reconstruction  based  on  numerically  generated  frequency  domain  data  with  0.1%  noise  in 
amplitude  and  0. 1  degree  noise  in  phase  is  demon strated.  The  size  of  the  domain  was  20cmxl0cmx4cm  and  the  grid 
was  81x41x33,  equivalent  to  0.25cmx0.25cmx0. 125cm  resolution.  The  size  of  the  inverse  grid  was  19x9x4, 
equivalent  to  1cm  voxels  and  1368  unknowns.  There  are  five  sources  with  the  same  frequency  modulation 
(140MHz)  and  81  detectors,  equivalent  to  810  measurements.  The  source  locations  are  (0,0,0),  (-1,1,0),  (1,1,0),  (-1,- 
1,0),  and  (1,-1, 0).  The  detector  locations  are  a  grid  from  -2  to  2  in  x  and  y  every  0.5cm.  There  are  two  lcm-cube 
heterogeneities  located  at  (0,1, 1.5)  and  (1,-1, 2.5).  The  optical  properties  are  pa=0.05cm‘1  and  Ps^-Ocm'1  for  the 
background,  p^O.lOcm'1  and  [is'=6.0cmx  for  the  cube  at  (0,1,1. 5),  and  (aa=0.08cm_1  and  fis’=10.0cm‘ 1  for  the  cube 
at  (1,-1, 2.5).  The  source  strength,  although  it  is  known  here,  is  generally  unknown  in  clinical  settings.  Therefore,  in 
order  to  remove  the  source  term  dependence  background  measurements  are  made  at  the  detector  locations  (<X>M°);  the 
corresponding  calculated  values  are  the  first  forward  solutions  (Oc°).  These  values  are  then  divided  out  of  the 
inverse  problem  for  each  measurement  (ie.  0M=0M/d>M°,  ^c^c^c)-  The  3-D  reconstruction  is  shown  in  Fig.  2 
along  with  the  actual  values.  Convergence  was  achieved  when  %2  changed  less  then  0.1%,  this  occurred  by  the  67th 
iteration.  The  reconstructed  cube  at  (0, 1,1.5)  had  pa=0. 08206cm*1  and  p,s*  =5. 9802cm'1  and  the  cube  at  (1,-1, 2.5)  had 
]Lia==0. 0772 lcm'1  and  ps,=12.2638cm'1.  The  total  time  was  5.6hrs  on  20nodes  (i.e.  1  master  and  19  slaves),  while 
synchronously  it  took  20.8hrs. 
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